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ABSTRACT 


This  report  summarizes  research  on  the  design  and  fatigue  analysis  of  laser-stake  welded 
connections  performed  at  the  University  of  Maine  from  January  2006  to  December  2007 
for  the  Structural  Response  of  Hybrid  Ship  Connections  Subjected  to  Fatigue  Loads 
(HSCF)  project.  The  HSCF  project,  funded  by  the  Office  of  Naval  Research,  is  part  of  a 
joint  effort  between  the  University  of  Maine,  Lehigh  University  of  Bethlehem,  PA  and 
Applied  Thermal  Sciences  Inc.  (ATS)  of  Sanford,  ME  an  engineering  research  and 
development  consulting  company  with  extensive  experience  in  implementation  of  laser 
welding  for  ships. 

Laser  welding  for  ship  structural  connections  is  a  maturing  technology  that  offers  much 
promise  for  increasing  productivity  and  reducing  fabrication  costs.  Past  fabrication 
methods  for  ship  structures  have  relied  upon  fastening  methods  such  as  periodically 
spaced  bolts  or  rivets,  adhesive  bonding  or  conventional  welding.  Laser  welding 
provides  an  extremely  robust  method  to  connect  steel  components  and  it  can  be 
performed  at  speeds  5-10  times  that  of  conventional  welding,  thus  requiring  less 
manufacturing  time.  Furthermore,  dimensional  accuracies  far  superior  to  that  obtained 
by  conventional  welding  have  been  realized  through  laser  welding.  Laser  welding  can 
also  be  used  in  metal/composite  hybrid  connections  wherein  modular  composite  panels 
with  embedded  steel  inserts  are  welded  to  a  steel  ship  hull. 

Since  ships  are  complex  structures  which  are  subjected  to  dynamic  loads,  a  thorough 
study  of  the  static  and  dynamic  load  resistance  is  essential  to  the  proper  design  of  laser 
welded  connections.  Accordingly,  the  major  objective  of  our  research  is  to  analyze  the 
static  and  fatigue  response  of  laser  welded  connections.  We  focus  on  continuous  laser 
welded  lap  joints  under  different  loading  conditions.  The  analyses  are  based  on  the 
principles  of  linear  elastic  fracture  mechanics  due  to  the  singular  stress  fields  at  the  edges 
of  the  laser  welds.  The  computational  model  presented  here  will  enable  engineers  to 
design  robust  laser  welded  lap  joints  for  prescribed  loads. 
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A  static  stress  analysis  is  performed  to  obtain  the  stress  intensity  factors  and  the  failure 
load  of  laser  welded  lap  joints.  The  geometry  and  nonlinear  contact  conditions  of  lap 
joints  make  it  difficult  to  obtain  an  analytical  closed-form  solution.  Therefore,  the 
commercial  finite  element  software  ABAQUS  is  used  to  perform  a  stress  analysis  of 
various  joint  configurations  and  obtain  the  stress  intensity  factors  and  joint  strength.  It  is 
found  that  the  overlap  length,  weld  width  and  plate  thicknesses  play  a  major  role  in 
determining  the  joints  behavior  under  different  loading  conditions.  The  results  are 
utilized  to  construct  a  metamodel  that  can  be  used  to  estimate  the  fracture  strength  of  lap 
joints  of  arbitrary  dimensions  without  having  to  perform  a  detailed  finite  element  analysis 
every  time.  The  parameters  in  the  metamodel  are  optimized  using  a  real-coded  single 
objective  genetic  algorithm.  The  metamodels,  which  are  validated  and  found  to  be  quite 
accurate,  are  very  useful  for  designing  lap  joints  for  prescribed  static  loads. 

The  fatigue  response  of  laser  welded  lap  joints  is  also  investigated  for  different  weld 
geometries  and  loading  conditions.  Fatigue  failure  is  often  a  complex  process  due  to 
numerous  factors  that  influence  crack  growth.  Nevertheless,  it  is  important  to  develop  a 
fatigue  model  for  purposes  of  scheduling  inspection  at  regular  intervals,  determining 
required  maintenance,  conducting  a  failure  analysis,  and  minimizing  the  problem  of  crack 
propagation  in  future  designs.  A  fracture  specific  finite  element  program  FRANC2D  is 
used  to  analyze  crack  growth,  obtain  the  crack  propagation  trajectories  and  compute  the 
corresponding  stress  intensity  factor  histories  for  the  lap  joint  geometries  under 
consideration.  The  Paris  fatigue  growth  law  is  used  to  obtain  the  S-N  curves  for  constant 
amplitude  cyclic  loads.  It  is  found  that  the  weld  width  has  a  significant  effect  on  the 
crack  trajectory  and  fatigue  life. 
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1.  INTRODUCTION 

This  report  for  the  Structural  Response  of  Hybrid  Ship  Connections  Subjected  to 
Fatigue  Loads  (HSCF)  project  covers  research  on  the  design  and  fatigue  analysis  of  laser- 
stake  welded  connections  performed  at  the  University  of  Maine  from  January  2006  to 
December  2007.  In  this  effort,  the  University  of  Maine  (UMaine)  is  partnered  with 
Lehigh  University  of  Bethlehem,  PA  and  Applied  Thermal  Sciences  Inc.  (ATS)  of 
Sanford,  ME  an  engineering  research  and  development  consulting  company  with 
extensive  experience  in  implementation  of  laser  welding  for  ships. 

The  U.S.  Navy  currently  has  an  objective  to  investigate  the  use  of  lightweight 
materials  to  enhance  the  future  naval  capabilities.  Composite/metal  hybrid  systems  offer 
a  solution  that  meets  this  criterion.  The  HSCF  program  at  the  University  of  Maine  is 
focused  on  the  study  of  the  fatigue  response  of  combined  composite  and  metal  systems 
and  innovative  connection  methods  such  as  laser  welding.  The  ultimate  goal  is  to 
achieve  an  efficient,  cost-effective  durable  ship  hull.  An  efficient  hull-form  structure 
must  be  lightweight  and  stiff  if  it  is  to  maintain  its  shape  while  resisting  the  applied  loads. 
It  also  must  be  fatigue,  impact  and  shock  resistant.  One  of  the  primary  cost  drivers  in 
developing  advanced  hull-forms  with  conventional  techniques  is  in  the  metal  forming  of 
complex  shapes.  Composite  materials  have  an  inherent  ability  to  perform  complex 
shaping  at  relatively  little  incremental  cost  compared  to  flat  panels.  Material  cost  is  a 
strong  consideration  in  Navy  ships,  which  are  massive  complex  structures  with  large 
amounts  of  material  used  in  their  hulls.  For  this  reason,  when  choosing  composites  for 
ships,  the  relatively  inexpensive  fiberglass  reinforced  polymer  (GRP)  composite  systems 
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using  a  vacuum  assisted  resin  transfer  molding  (VARTM)  process  are  feasible  and  much 
research  effort  has  focused  on  this  material  system. 

The  long-term  goal  of  the  research  performed  at  UMaine  is  to  develop  and 
implement  robust  techniques  for  evaluation  of  the  fatigue  and  environmental  response  of 
hybrid  construction  for  use  in  U.S.  Navy  ships.  Current  design  goals  for  the  new  class  of 
navy  ships  are  for  vessels  that  are  much  faster  and  stealthier  than  currently  available  in 
the  fleet.  These  newer  vessels  will  dictate  the  use  of  advanced  materials  in  their  design 
that,  among  other  characteristics,  are  lightweight  and  nonmagnetic.  Although  the  original 
direction  for  use  of  new  materials  focused  on  composites,  it  was  quickly  determined  that 
glass  reinforced  polymer  (GRP)  composites  alone  lack  the  stiffness  and  connectability 
required  in  the  medium  to  long  length  modem  ships.  In  these  circumstances  the  use  of  a 
hybrid  composite/metal  structure  can  combine  the  benefits  of  traditional  steel 
construction  while  using  the  low-density  nature  of  composites  in  weight  critical  areas. 
Recently,  the  application  of  hybrid  metal-to-composite  structures  is  found  in  a  widening 
number  of  engineering  disciplines  due  to  the  inherent  benefits  of  such  systems. 

A  ship  is  a  dynamically  loaded  vessel  and  mitigation  of  structural  failures  is 
essential.  More  often  than  not  structural  failures  occur  due  to  fatigue  at  connections  and 
interfaces.  Therefore,  study  of  connection  fatigue  is  essential  in  a  robust  structural  design. 
An  accurate  appraisal  of  structural  integrity  depends  primarily  on  proper  assessment  of 
the  structural  response  of  the  connections  and  interfaces  and  a  sound  estimate  of  the 
dynamic  nature  of  the  loads  that  induce  failure.  Accordingly,  one  must  perform  a 
thorough  investigation  into  the  mechanics  of  the  connections  and  interfaces  of  the  vessel. 
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Understanding  connection  response  and  development  of  connections  with  highly  efficient 
structural  response  is  an  essential  part  of  this  research. 

Shipbuilders  are  showing  increasing  interest  in  laser  welding  as  a  means  of 
reducing  fabrication  costs  and  distortion  in  ship  construction  and  repair.  Current  interest 
lies  in  laser  welded  fabricated  shapes  and  Laser  welded  Sandwich  panels.  It  is  ultimately 
envisioned  that  laser  welding  will  be  as  common  in  the  shipyard  as  other  processes  such 
as  MIG,  TIG  and  SMAW.  Laser  stake-welding  of  the  outer  shell  of  a  ship  structure  to  its 
frame  is  depicted  in  Fig.  1  (a).  Laser  welding  can  also  be  used  in  metal/composite  hybrid 
connections.  In  this  method,  modular  composite  panels  with  embedded  steel  inserts  are 
welded  to  a  steel  ship  hull  as  shown  in  Fig.  1(b).  However,  the  high  heat  input  from 
conventional  welding  techniques  can  be  detrimental  to  the  polymer  matrix  composite 
material.  In  comparison,  the  laser  welding  process  allows  welding  to  occur  much  closer 
to  the  composite  material  due  to  localized  heating. 

Cao  and  Grenestedt  (2003)  and  Clifford  et  al.  (2002)  performed  a  study  on  a 
particular  hybrid  joint  ideal  for  laser  welding.  In  order  to  improve  the  strength  and 
mechanical  interlock  properties  of  embedded  hybrid  joints,  Melograna  and  Grenestedt 
(2002)  investigated  the  effect  of  perforations  in  the  steel  on  the  strength  of  such  joints. 
Joints  were  manufactured  in  which  steel  inserts  with  a  hole  pattern  patented  by  Unden 
and  Ridder  (1985)  were  laminated  into  the  GRP.  It  was  found  that  the  use  of  perforations 
significantly  aided  the  strength  of  the  joints.  There  were  multiple  reasons  for  this, 
including  reduction  in  stiffness  mismatch  in  the  elastic  properties  of  the  steel  and 
composite  attributed  to  the  progressively  widening  hole  pattern  in  the  steel.  Additionally, 
the  holes  provided  significant  mechanical  interlock  whereas  embedded  joints  without 
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holes  rely  strictly  on  the  adhesion  between  the  insert  and  the  composite  for  strength. 
Understanding  the  dynamic  response  of  laser  welded  connections  in  hybrid 
composite/metal  connections  is  critical  if  they  are  to  be  used  successfully  in  ship 
building. 

The  major  objective  of  the  present  research  is  to  evaluate  the  static  and  fatigue 
response  of  laser  welded  connections.  We  focus  on  continuous  laser  welded  lap  joints 
under  different  loading  conditions.  A  static  stress  analysis  is  performed  first  to  obtain  the 
stress  intensity  factors  and  the  failure  load  of  laser  welded  lap  joints.  Subsequently,  the 
fatigue  response  of  laser  welded  connections  is  studied  using  a  fracture  specific  finite 
element  model.  The  computational  model  presented  here  will  enable  engineers  to  design 
robust  laser  welded  lap  joints  for  ship  structures. 
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(a) 


(b) 


Figure  1.1:  Laser  stake  welds  in  (a)  ship  construction  (b)  metal/composite  hybrid 

connection. 
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1.1  Background  and  Research  Motivation 

In  recent  years,  high  power  density  welding  technologies  like  laser  welding  have 
been  increasingly  utilized  in  industrial  manufacturing.  This  is  due  to  localized  energy 
input,  increased  welding  speed,  high  penetration,  easy  automation,  excellent  accuracy  and 
small  heat  affected  zones.  These  qualities  lead  to  smaller  dimensional  and  shape 
distortion  of  the  welded  pieces.  With  this  in  mind,  the  task  of  analyzing  laser  stake 
welded  connections  to  determine  their  load  carrying  capacity  and  fatigue  response  is 
undertaken.  In  addition,  the  information  obtained  by  performing  various  failure  analyses 
is  cleverly  utilized  to  design  optimal  laser  welded  joint  configurations. 

1.1.1  Laser  Welding 

The  term  LASER  is  an  acronym  for  "Light  Amplification  by  Stimulated  Emission 
of  Radiation"  and  is  defined  as  a  device  that  emits  highly  amplified  and  coherent 
radiation  of  one  or  more  discrete  frequencies.  In  the  laser  welding  process,  the  photons 
are  targeted  on  the  work  piece  surface  which  needs  to  be  welded.  Highly  concentrated 
photon  energy  gets  converted  into  thermal  energy,  which  causes  the  materials  to  melt  and 
fuse  together  as  the  laser  energy  is  removed.  Laser  beam  welding  (LBW)  is  a  type  of 
fusion  welding  that  is  used  to  produce  selective  area  spot  welds  or  linear  continuous  seam 
welds.  LBW  is  capable  of  welding  carbon  steels,  HSLA  steels,  stainless  steel,  aluminum, 
and  titanium.  Due  to  high  cooling  rates,  cracking  is  a  concern  when  welding  high-carbon 
steels.  The  weld  quality  is  high,  and  is  similar  to  that  of  electron  beam  welding.  High 
power  gas  lasers  are  especially  suitable  for  high  volume  applications,  and  as  a  result, 
LBW  is  rather  prevalent  in  the  automotive  industry. 
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There  are  two  types  of  laser  welding  processes,  conduction  and  penetration. 
Laser  conduction  welding  relies  on  the  conductivity  of  a  material.  It  is  used  for  spot, 
continuous  and  partial  penetration  seam  welding.  Laser  penetration  welding,  on  the  other 
hand,  is  carried  out  by  focusing  the  laser  beam  energy  at  a  single  location  until  the 
abutted  materials  are  heated  to  a  liquid  state.  The  laser  energy  is  removed  and  the  liquid 
material  begins  to  cool,  thus  joining  the  previously  separate  pieces  as  the  liquid  metal 
solidifies.  Another  important  process  is  hybrid  laser  arc  welding  (HLAW)  which 
combines  laser  with  conventional  gas  metal  arc  welding  processes.  This  increases  the 
speed  of  performing  the  weld  and  provides  higher  strength  bonds  with  less  potential  for 
undercutting. 

The  two  types  of  lasers  commonly  used  in  metal  working  are  solid-state  lasers 
and  gas  lasers  (especially  CO2  lasers).  The  first  uses  any  solid  media,  including  synthetic 
ruby  and  chromium  in  aluminum  oxide,  neodymium  in  glass,  and  the  most  common  type 
being  crystal  YAG  (composed  of  yttrium,  aluminum,  and  garnet)  doped  with 
neodymium.  Gas  lasers  use  mixtures  of  gases  like  helium,  nitrogen,  and  carbon  dioxide 
as  a  medium.  Solid-state  lasers  operate  at  wavelengths  of  1 .06  pm,  much  shorter  than  gas 
lasers.  Both  lasers  can  operate  in  both  pulsed  and  continuous  modes.  Typical  power 
output  for  ruby  lasers  is  10-20  W  while  the  YAG  laser  outputs  between  0.04  and  600  W. 
To  deliver  the  laser  beam  to  the  weld  area,  fiber  optics  is  usually  employed.  Gas  lasers 
pertaining  to  this  project,  use  high-voltage,  low-current  power  sources  to  supply  the 
energy  needed  to  excite  the  gas  mixture  used  as  a  lasing  medium.  These  lasers  have 
wavelengths  of  10.6  pm.  As  there  are  no  waveguides  used  for  gas  systems,  a  lens  and 
mirror  delivery  system  is  used.  Power  outputs  for  gas  lasers  can  be  much  higher  than 
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solid-state  lasers,  reaching  25  kW.  It  is  hopeful  that  modular  fiber  lasers  can  achieve  the 
same  power  output. 

1.1.2  Stresses  due  to  Welding 

The  technique  of  laser  welding,  which  can  join  steel  plates  at  speeds  of  several 
meters  per  minute  in  a  single  weld  pass,  is  very  promising.  This  type  of  technology  could 
prove  very  useful  for  huge  structures  like  naval  ships,  oil  tankers,  etc.  However,  the  fast 
cooling  rates  in  the  welded  zone  produce  a  very  hard  and  brittle  microstructure  which 
may  produce  structural  integrity  problems.  This  area  is  called  the  heat-affected  zone 
(HAZ).  The  size  and  strength  of  the  HAZ  depends  on  the  material’s  thermal  properties 
and  the  power  and  speed  of  weld  process  being  used.  The  thermal  diffusivity  of  the  base 
material  plays  a  large  role.  If  the  diffusivity  is  high,  the  material  cooling  rate  is  high  and 
the  HAZ  is  relatively  small.  Conversely,  a  low  thermal  diffusivity  leads  to  slower  cooling 
and  a  larger  HAZ.  Laser  beam  welding  gives  a  highly  concentrated  and  limited  amount  of 
heat,  resulting  in  a  small  HAZ,  which  is  desirable.  While  welding,  the  heat  flows  from 
high  temperature  metal  to  the  surrounding  metal  which  generates  thermal  gradients 
causing  thermal  stresses.  When  this  material  cools  and  solidifies  it  retains  residual 
stresses.  There  are  processes  available  to  alleviate  these  residual  stresses,  for  example, 
annealing.  In  this  work,  any  connection  being  analyzed  is  assumed  to  have  been  already 
stress  relieved. 

1.1.3  Fracture  Mechanics 

Fracture  mechanics  is  comprised  of  a  series  of  models  used  to  describe  the 
influence  of  cracks  or  defects  on  material  failure.  The  aim  for  a  portion  of  this  work  is  to 
use  fracture  mechanics  to  analyze  a  laser  welded  lap  joint  model  under  static  and  fatigue 
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loads.  To  demonstrate  the  need  for  doing  this,  the  joint  is  modeled  using  a  commercial 
finite  element  package  ABAQUS.  When  certain  unit  moments  are  applied  at  the  free  ends 
of  the  welded  plates,  it  is  observed  that  the  joint  opens  up  and  stresses  at  the  weld-ends 
start  approaching  infinity.  Thus,  it  is  reasonable  to  say  that  the  non-welded  overlapped 
lengths  act  as  cracks  or  notches  and  the  weld  ends  as  left  and  right  crack  tips. 

Griffith  (1920)  performed  fundamental  experiments  on  the  strength  of  thin  glass 
fibers.  He  concluded  that  thinner  fibers  failed  at  higher  stresses  than  thicker  fibers,  an 
outcome  not  predicted  by  conventional  stress-based  strength  theory.  Griffith  considered 
the  thermodynamic  equilibrium  of  a  crack  system  and  stated  that  a  crack  will  grow  when 
the  total  energy  of  the  system  starts  reducing.  In  a  classic  paper,  Inglis  (1913),  came  up 
with  an  expression  in  which  the  magnitude  of  stress  at  the  tip  of  a  hole  varies  both  with 
the  size  of  the  hole  and  sharpness  of  the  tip,  thus  providing  a  first  means  in  determining 
the  larger  stress  gradients  around  the  crack  tips. 

Irwin  (1956)  developed  an  energy  release  rate  concept  which  is  an  extension  of 
the  Griffith’s  theory,  albeit  in  a  more  useful  form.  Before  that,  Westergaard  (1939) 
developed  a  semi-inverse  technique  for  analyzing  stresses  and  displacements  ahead  of  a 
sharp  crack.  Irwin  used  Westergaard’ s  approach  to  show  that  the  stresses  and 
displacements  near  the  crack  tip  could  be  described  by  a  single  parameter  that  was  related 
to  the  energy  released  rate,  which  is  called  as  the  stress  intensity  factor  (SIF).  Over  time, 
there  have  arisen  two  fundamental  approaches  to  fracture  problems;  the  strain  energy 
release  rate  method,  which  is  based  on  global  energy  balance  of  the  whole  system,  and 
the  stress  intensity  factor  method,  which  is  based  on  the  local  stress  distribution  around  a 
crack  tip.  For  the  sake  of  ease,  simplicity,  clear  demarcation  of  the  various  fracture 
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modes  and  their  individual  effects  on  the  failure  strength  of  the  joint,  and  ultimately  to  be 
able  to  develop  a  simple  method  of  predicting  fatigue  behavior  of  the  joints,  the  stress 
intensity  factor  approach  will  be  used  throughout  this  work. 


1.1. 3.1  Stress  Intensity  Factor 

The  traditional  problem  with  the  strain  energy  release  rate  method  is  that  it 
depends  on  the  total  energy  state  of  the  whole  structure  before  and  after  the  crack  growth. 
For  complex  problems  which  typically  involve  large  models  and  numerical  solution 
techniques,  this  approach  is  not  very  feasible.  As  a  result,  the  stress  intensity  factor 
method  can  be  used  as  a  simple  to  implement  alternative. 

Consider  the  tip  of  a  crack  in  a  uniform  tension  field  o,  as  shown  in  Fig.  1 .2.  In  a 
polar  coordinate  system,  the  general  two  dimensional  solution  for  the  stresses  near  the  tip 
of  the  crack  at  an  angle  6  and  radial  distance  r  are: 

K,  0  .  e  .  30. 

<x  =  .  cos— (1  -  sin  — sin — ), 

■N &  2 V  2  2  7 

k.  e  n  .  e  .  30.  (i) 

cr  =  . — . —  cos — (1  +  sin— sin — ), 

'  yflm-  2 V  2  2  ' 

K.  .  0  0  30 

r  =  --J=zz  sin  —  cos  —  cos  — . 

*  ^  2  2  2 

Note  that  each  stress  component  is  proportional  to  Ki  which  is  known  as  mode-I  stress 
intensity  factor.  K/  can  be  defined  as  the  intensity  of  the  stress  singularity  that  occurs  at 
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Figure  1.2:  Crack  tip  in  a  uniform  tension  field. 
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the  tip  of  a  sharp  crack  due  to  a  pure  transverse  loading.  If  Kj  is  known,  then  the  entire 
stress  distribution  at  the  crack  tip  can  be  computed  by  equation  (1).  Kj  can  thus  be  said  is 
a  function  of  the  applied  load,  the  crack  length  and  the  specimen’s  geometry.  Using  the 
concept  of  dimensional  analysis,  K/  can  be  evaluated  as: 

K,=P<r4d,  (2) 

where  ft  is  a  dimensionless  geometry  parameter,  c  is  the  far  field  stress  magnitude  and  cl 
is  the  crack  length.  The  unit  of  stress  intensity  factor  is  PsiVin.  Because  K )  is  defined  as  a 
measure  of  stress  intensity  at  the  crack  tip,  it  can  be  said  that  a  crack  will  grow  when  K/ 
reaches  a  critical  value  K )  =  K/c,  where  K/c  is  the  critical  stress  intensity  factor  for  a 
material.  K/c  is  thus  considered  a  material  property  (analogous  to  yield  strength)  that 
defines  the  material’s  resistance  to  crack  growth.  It  is  often  referred  to  as  fracture 
toughness  of  the  material. 

Three  modes  of  fractures  can  be  defined;  the  previously  mentioned  mode-I  which 
is  due  to  in-plane  tension,  mode-II  due  to  in-plane  shear  and  mode-III  where  there  is  out 
of  plane  shear  load.  These  modes  are  illustrated  in  Fig.  1.3.  Mode-I  is  typically  seen  as 
the  dominant  mode  in  many  structural  materials  and  receive  more  attention  in  literature. 
Each  mode  has  its  own  independent  critical  stress  intensity  factor,  these  being  K/c,  Kl/c 
and  Kmc  respectively. 

1.1. 3.2  Other  Fracture  Criteria 

The  basic  assumption  of  Linear  Elastic  Fracture  Mechanics  (LEFM)  is  that  all 
available  strain  energy  goes  into  propagating  a  crack  or,  according  to  Griffith,  in  creating 


12 


Mode  II  -  In-plane  shear 


Figure  1.3:  Illustration  of  the  three  fracture  modes. 
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new  surfaces.  In  reality,  however,  many  materials  have  some  micro  structural 
mechanisms  that  dissipate  the  available  energy,  for  instance: 

•  Plastic  deformation  around  the  crack  tip  in  metals. 

•  Micro-cracking  and  friction  in  concrete  and  rock. 

•  Fiber  bridging  in  wood  and  fibrous  composites. 

•  Inclusions. 

These  phenomena  can  be  referred  to  as  toughening  mechanisms.  The  only 
toughening  mechanism  seen  in  isotropic  metals  is  plastic  deformation.  Although  in 
metals,  LEFM  concepts  can  be  slightly  altered  in  order  to  tolerate  limited  plasticity  in  the 
crack  tip  region,  the  validity  of  LEFM  ceases  when  significant  plastic  deformation 
occurs.  At  high  loads  in  a  lap  joint  configuration  made  of  steel,  it  is  safe  to  assume  that 
the  energy  loss  is  due  to  plastic  deformation.  Beyond  the  elastic  limit,  another  concept, 
Elastic-Plastic  Fracture  Mechanics  (EPFM)  comes  into  the  picture.  Although  it  is  an 
extension  to  the  elastic  deformation,  EPFM  cannot  be  applied  for  yield  leading  to  plastic 
collapse.  Wells  (1961)  attempted  to  apply  LEFM  to  low  and  medium  carbon  steels  and 
noticed  the  displacement  of  crack  surfaces  near  the  crack  tip  due  to  plastic  deformation 
preceding  failure.  As  a  result,  he  proposed  an  alternative  fracture  criterion  known  as 
Crack  Tip  Opening  Displacement  (CTOD). 

Rice  (1968)  in  a  seminal  paper  idealized  plastic  deformation  as  a  nonlinear  elastic 
process  and  developed  another  parameter  to  characterize  nonlinear  material  behavior 
ahead  of  a  crack.  He  expressed  energy  release  rate  for  non-linear  materials  as  a  line 
integral  which  he  termed  as  the  J-Integral.  The  J-Integral  can  be  evaluated  along  any 
arbitrary  contour  around  the  crack  tip  as  it  is  path  independent. 
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Linear  elastic  fracture  mechanics  has  proven  to  be  very  useful  for  the  fatigue 
process  in  the  near-tip  region  of  a  crack  under  cyclic  loading.  The  idea  of  using  LEFM  is 
that  the  stress  intensity  factor  can  still  be  used  to  analyze  progressive  crack  propagation 
in  the  fatigue  process  zone.  In  the  case  of  a  constant  amplitude  cyclic  load,  experimental 
studies  have  shown  that  the  crack  growth  rate  is  a  function  of  the  amplitude  of  the  stress 
intensity  factor.  One  such  fatigue  crack  growth  law  was  proposed  by  Paris  and  Erdogan 
(1960).  The  theoretical  predictions  of  the  Paris  Law  compare  favorably  with  experimental 
results  for  moderate  loads  considered  in  this  research.  Realistically,  micro  structural 
mechanisms  (like  plasticity)  also  affect  the  prediction  of  the  fatigue  life  of  a  joint.  This 
requires  the  consideration  of  an  elasto-plastic  failure  criterion  into  the  fatigue  equation. 
Since  such  an  analysis  is  complicated,  the  present  research  utilizes  a  simpler  and 
reasonably  correct  approach  based  on  LEFM  to  outline  the  procedure  of  predicting 
fatigue  life  of  elastic  and  moderately  plastic  materials. 

1.1.4  Single-objective  Genetic  Algorithms 

The  mathematical  model  for  the  computation  of  the  fracture  strength  of  arbitrary 
joint  dimensions  is  optimized  using  real-coded  single  objective  GA.  Several  techniques 
have  been  created  for  optimizing  a  given  objective  function. 
These  include  the  linear  least  squares  method,  the  steepest  descent  method,  the 
Levenberg-Marquardt  method,  and  the  Gauss-Newton  method  among  many  others.  These 
methods,  unlike  the  genetic  algorithm  (GA),  tend  to  hone  in  on  local  minima  or  maxima 
and  are  only  efficient  when  the  objective  function  possesses  continuous  derivatives  and  is 
not  too  complex.  Genetic  algorithms  are  used  when  the  search  space  is  either  too  large  or 
too  complex  for  an  analytical  method.  Genetic  algorithm,  an  optimization  methodology. 
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introduced  and  investigated  by  Holland  (1975),  is  inspired  by  evolutionary  processes  and 
the  “survival  of  the  fittest”  concept.  In  broader  terms,  a  GA  is  a  stochastic  global  search 
method  that  mimics  the  concept  of  natural  biological  evolution. 

To  begin  the  GA  procedure,  a  solution  space  is  first  provided  and  a  population  of 
random  potential  solutions  is  generated.  The  fitness  of  each  individual  is  evaluated  and  a 
fitness  rating  is  assigned  to  individuals  reflecting  the  designer’s  goals.  Applying  the 
principle  of  “survival  of  the  fittest”,  this  fitness  rating  is  used  to  select  and  keep  multiple 
copies  of  better  performing  individuals  and  eliminate  poor  design  solutions.  Individuals 
are  chosen  from  among  the  better  solutions  retained,  to  “mate”  and  create  child  solutions 
using  recombination  operators  borrowed  from  natural  genetics,  namely  crossover  and 
mutation,  which  in  general  leads  to  a  child  population  of  mostly  improved  designs.  The 
new  child  population  is  often  compared  with  the  previous  parent  population  and  only  the 
best  of  the  parent  and  child  populations  is  retained.  This  process  is  called  elitism  and  is 
often  used  for  increasing  the  performance  of  the  genetic  algorithm.  The  process  of 
selecting  the  best  individuals  for  mating,  creating  child  offspring  populations  and  so  on  is 
repeated  over  many  generations  until  one  is  satisfied  with  the  final  result  or  some 
established  convergence  criteria  is  met.  Fig.  1.4  illustrates  the  GA  process  in  a  simple 
flowchart. 
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Figure  1.4:  Genetic  algorithm  procedure. 
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1.2  Literature  Review 


There  is  an  increasing  interest  by  ship  builders  to  incorporate  laser  welding 
technology  in  their  operations  to  reduce  fabrication  costs.  Before  one  can  start  using  laser 
welding  processes  in  the  construction  of  ship  structures,  one  must  first  know  the 
durability  or  fracture  toughness  and  fatigue  resistance  of  a  laser  welded  structure. 
Sumpter  (1996),  studied  the  cause  of  fractures  in  ship  structures  on  the  microscopic  level. 
The  quality  of  laser  welding  is  discussed  in  the  work  which  gamers  its  insights  with  the 
help  of  experiments  that  place  cracks  at  different  locations  in  the  structure.  These 
locations  include  cracks  at  the  weld  centerline,  at  heat  affected  zone  or  in  the  plate 
material.  The  obtained  results  are  then  compared  with  those  obtained  by  similar 
experiments  where  arc  welding  was  used.  During  the  tests  the  major  objective  was  to 
avoid  cleavage  fracture,  the  occurrence  of  which  is  strongly  influenced  by  temperature. 
The  toughness  of  a  laser  welded  mild  steel  ship  plate  was  measured  by  dynamic  fracture 
mechanics  testing  at  a  temperature  of  0°C.  The  elastic-plastic  fracture  criterion  is  found 
to  be  the  smallest  when  the  crack  was  located  at  the  HAZ,  and  is  similar  in  numerical 
value  as  compared  to  material  joined  by  typical  arc  welding.  All  in  all,  the  work  provided 
helpful  insights  on  the  process  of  determining  K/c  of  laser  welded  materials. 

Caccese  et  al.  (2006)  compared  the  effect  of  weld  geometries  on  laser  welded 
HSLA-65  steel  cruciforms.  New  welding  techniques  were  employed,  specifically  a 
combined  laser  and  gas-metal-arc  welding  (GMAW)  procedure.  A  number  of  fatigue 
tests  were  run  on  laser- welded  cruciforms  of  different  weld  geometries  available  with 
laser  welding,  as  well  as  some  conventionally  welded  specimens.  It  was  found  that  if 
weld  geometry  is  carefully  controlled,  the  hybrid  laser  welded  specimens  performed 
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better  than  the  conventionally  welded  specimens  and  what  the  currently  accepted  fatigue 
life  equations  would  predict. 

In  a  typical  lap  joint  geometry,  it  can  be  said  that  there  exists  two  interface  cracks 
between  the  two  plates.  An  approach  at  tackling  the  problem  is  given  by  Suo  and 
Hutchinson  (1990).  They  considered  a  semi-infinite  interface  crack  between  two  infinite 
isotropic  elastic  layers  under  general  edge  loading  and  presented  a  mathematical  model 
for  its  analysis.  Their  study  is  concentrated  mainly  on  interface  cracks  between  thin  films 
and  substrates.  They  have  introduced  expressions  for  the  stress  intensity  factors  for  the 
geometry  similar  to  the  geometry  and  loading  conditions  under  consideration. 
Furthermore,  their  computation  of  stress  intensity  factors  involved  the  principle  of 
superposition  and  thus,  the  technique  is  not  suitable  for  geometrically  non-linear 
structures  such  as  the  ones  we  consider  here  where  interface  contact  is  present.  Further 
discussion  of  this  topic  is  given  in  Section  2.6.4. 

When  the  crack  length  is  comparable  to  the  characteristic  dimensions  of  the 
structure,  it  becomes  difficult  to  compute  an  analytical  expression  for  the  stress  intensity 
factor  (Sun,  2001).  In  such  cases  numerical  techniques  become  very  useful.  One  such 
technique,  created  by  Irwin  (1956),  is  the  Crack  Closure  Integral  which  utilizes  nodal 
force  and  displacement  information  at  the  crack  tip  to  compute  the  stress  intensity  factors. 
This  force  and  displacement  data  can  now  easily  be  obtained  from  any  finite  element 
package.  Several  researchers  have  found  the  technique  useful,  but  a  little  onerous.  The 
method  requires  two  separate  analyses  to  obtain  the  desired  results,  one  with  the  original 
geometry  to  get  the  nodal  forces  and  a  second  with  an  extended  crack  to  get  the 
displacements.  Improving  upon  this  concept,  Rybicki  and  Kanninen  (1977)  came  up  with 
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a  more  efficient  technique  for  evaluating  stress  intensity  factors  using  finite  element 
calculations  called  the  Modified  Crack  Closure  Technique.  The  method  falls  under  the 
category  of  indirect  methods  and  permits  evaluation  of  both  mode-I  and  II  stress  intensity 
factors  from  the  results  of  a  single  analysis.  In  their  work  they  define  the  limitations  of 
the  method  and  a  correction  factor  for  the  case  where  unequal  element  lengths  are  present 
before  and  after  the  crack  tip.  This  method  proved  to  be  very  helpful  in  this  work  with 
regard  to  computing  stress  intensity  factor  values  for  different  joint  dimensions. 

Unlike  Rybicki  and  Kanninen,  Krueger  (2002)  referred  to  the  aforementioned 
technique  as  the  Virtual  Crack  Closure  Technique  (VCCT)  and  presented  it  with 
reference  to  the  delamination  of  composite  structures  and  used  the  method  to  compute 
damage  tolerance.  In  addition,  Krueger  also  presented  several  ways  in  which  the  VCCT 
can  be  applied.  For  example,  Krueger  has  supplied  equations  for  two-dimensional 
quadrilateral  elements  with  linear  and  quadratic  shape  functions,  formulae  for  using  the 
VCCT  method  with  three-dimensional  solid  elements  and  also  expressions  for  using  the 
method  with  plate/shell  elements.  While  the  correction  factor  for  unequal  element  length 
was  already  established,  he  has  presented  an  analogous  correction  for  unequal  element 
widths. 

A  work  that  uses  a  very  similar  approach  to  the  one  adopted  here,  is  performed  by 
Gowhari  and  Anaraki  (2005).  They  analyzed  the  fatigue  response  of  spot-welded  lap 
joints  for  shear-tension  loading.  A  combination  of  EPFM  and  LEFM  techniques  were 
used  and  results  were  presented  in  the  form  of  structure,  notch  and  crack  configuration 
factors.  Their  main  aim  was  to  present  fracture  data  in  terms  of  non-dimensional 
geometric  parameters  in  order  to  provide  a  direct  method  for  weld  life  prediction.  In  their 
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work  they  also  assumed  the  non-welded  overlapped  plate  length  as  a  notch,  and  since  the 
welds  under  consideration  resembled  a  nugget,  the  notch  was  considered  a  three 
dimensional  crack.  A  finite  element  parametric  study  was  performed  which  used 
normalized  forms  of  some  of  the  influential  geometric  dimensional  parameters. 
Predictive  equations  for  the  SIF  and  the  J-integral  are  given  by  employing  a  statistical 
multiple  non-linear  regression  model  on  the  results  obtained  from  the  parametric  study. 
Unlike  our  work,  they  used  three-dimensional  finite  element  models  and  used  over  1000 
analyses  to  obtain  accurate  prediction  equations. 

Terasaki  and  Kitamura  (2004)  presented  a  procedure  for  prediction  of  static 
fracture  strength  of  laser  welded  lap  joints  by  numerical  analysis  is  presented.  In  their 
work,  the  finite  element  model  is  compared  with  the  experimental  results  at  several  steps, 
even  to  decide  the  appropriate  boundary  conditions  for  the  model.  They  have  established 
a  numerical  method  which  conforms  to  the  experimental  load-displacement  curves  of 
tensile-shear  tests,  and  then  used  that  method  to  investigate  whether  the  experimentally 
determined  static  fracture  strength  of  lap  joints  can  be  predicted  by  the  equivalent  plastic 
strain.  The  analyses  were  performed  on  ABAQUS  with  CPS8R  (8-node  two  dimensional 
reduced  integration)  elements  using  large  deformation  theory  where  plane-stress 
conditions  prevail. 

The  analytical  expressions  for  standard  cracked  geometries,  like  the  edge  crack 
and  center  crack  models  with  tensile  and  shear  loading,  are  helpful  for  assessing  the 
accuracy  of  finite  element  models  by  performing  convergence  studies  that  compare  the 
analytic  and  numerical  results.  The  analytical  expressions  for  stress  intensity  factors  for 
several  standard  configurations  are  provided  in  books  by  Broek  (1978)  and  Irwin  (1985). 
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A  very  clear  explanation  of  the  mixed  mode  plastic  zone  size  is  given  in  a  book  by 
Gdoutos  (2005),  both  for  plane-stress  and  plane-strain  conditions.  Crack  growth,  elastic- 
plastic  fracture  mechanics  and  fatigue  have  also  been  discussed  in  detail,  some  topics  of 
which  are  pursued  in  this  work  in  Chapter  3. 

Assessing  the  strength  of  a  connection  design  requires  not  only  analyzing  the 
joints  resistance  to  static  loads  but  also  for  determining  its  response  to  fatigue  conditions 
as  well.  Recent  experiments  in  the  fatigue  performance  of  laser  welded  lap  joints  show 
that  the  performance  of  the  connection  is  dependent  upon  the  weld  pattern.  Multiple 
longitudinal  welds  possess  higher  fatigue  strength  than  simple  transverse  welds.  Hsu  and 
Albright  (1991)  established  a  technique  for  predicting  the  fatigue  response  of  laser 
welded  lap  joints.  A  static  stress  analysis  was  combined  with  a  cyclic  stress-strain 
analysis  to  establish  a  model  for  calculating  joint  stresses  under  cyclic  loading. 

Laser  welded  lap  joints  suffer  from  defects  like  notch  effects,  surface  cracks, 
residual  stresses,  etc.  Cho  et  al.  (2004)  performed  finite  element  analyses  and  fatigue  tests 
to  estimate  the  fatigue  strength  of  laser  welded  lap  joints.  They  determined  that  residual 
stresses  and  welding  direction  were  two  major  factors  that  affected  the  fatigue  strength  of 
the  joint.  When  a  fatigue  load  is  applied  to  a  laser  welded  joint  containing  residual 
stresses  the  applied  and  residual  stresses  combine,  thus  reducing  the  fatigue  strength  of 
the  joint  significantly.  To  arrive  at  these  conclusions,  three  dimensional  models  with 
longitudinal  and  transverse  weld  directions  and  eight  node  quadrilateral  elements  are 
constructed  and  analyzed.  Using  temperature  distribution  data  from  a  heat  flow  analysis, 
a  thermo-elastic-plastic  stress  analysis  is  performed  which  simulates  the  complete  cooling 
process  all  the  way  down  to  the  ambient  temperature.  Next,  the  equivalent  fatigue 
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strength  is  calculated  using  the  elastic-plastic  stress  analysis  with  the  residual  stress  effect 
in  place.  The  work  was  performed  using  the  finite  element  program  ABAQUS  utilizing  a 
few  key  user-defined  subroutines. 

Computation  of  stress  intensity  factors  with  numerical  methods  is  typically 
carried  out  using  the  finite  element  methodology.  The  reader  is  referred  to  Logan  (2001) 
for  an  explanation  of  the  basic  and  applied  concepts  of  finite  element  analyses  which 
have  been  helpful  in  completing  this  work.  ABAQUS  documents  (Analysis  &  Theory 
manuals  Version  6.5)  have  been  extensively  referred  to  for  information  pertaining  to 
mesh  type,  mesh  control  features,  element  types  and  the  compatibility  between  the 
quadratic  and  triangular  elements. 

When  the  results  obtained  by  a  parametric  study  exhibit  a  strongly  non-linear 
behavior,  it  becomes  difficult  to  create  useful  predictions  by  using  simple  curve  fits  on 
the  data.  A  more  robust  technique,  namely  metamodels,  allows  one  to  construct  an 
approximate  model  using  function  values  at  a  few  sampling  points  typically  derived  from 
experimental  design  methods.  A  few  of  the  most  popular  metamodeling  techniques  are 
the  response  surface  methodology  (RSM),  the  kriging  method,  multivariate  adaptive 
regression  splines,  neural  networks  and  radial  basis  functions  (RBF).  Fang  et  al.  (2005) 
performed  a  comparative  study  of  two  metamodeling  methods,  RSM  and  RBF,  by 
considering  multi-objective  optimization  of  a  vehicle  body  subjected  to  a  frontal 
collision.  The  RSM  and  RBF  models  were  validated  via  finite  element  simulations  of  the 
full  scale  vehicle  model.  The  main  objective  of  the  work  is  to  assess  the  ability  of  the 
RSM  and  RBF  methods  as  applied  to  the  study  and  optimization  of  vehicle 
crashworthiness,  a  highly  non-linear  problem.  In  their  findings.  Fang  et  al.  note  that  RSM 
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may  not  be  suitable  for  modeling  highly  non-linear  responses  using  only  a  limited 
number  of  response  samples.  Radial  basis  functions,  on  the  other  hand,  were  found  to  be 
very  promising.  Although  the  RBF  models  are  computationally  more  expensive,  the 
optimization  results  were  found  to  be  more  accurate.  Based  on  the  findings  of  Fang  et  ah, 
we  utilize  the  augmented  radial  basis  function  method  to  construct  the  metamodel 
(Krishnamurth,  2003)  for  this  work.  Further  details  concerning  this  process  will  be 
discussed  in  Section  2.7. 

1.3  Outline  and  Research  Approach 

The  second  chapter  of  the  report  presents  a  method  for  the  evaluation  of  stress 
intensity  factors  at  the  notch  tips  in  laser  welded  lap  joints.  Throughout  the  work,  we 
limit  our  focus  to  steel  materials.  Residual  stresses  and  material  alterations  due  to  the 
laser  welding  process  have  been  neglected  in  the  present  analysis  for  the  sake  of 
simplicity.  We  assume  that  all  contacts  are  hard  and  ffictionless  and  confine  ourselves  to 
working  within  the  realm  of  linear  elastic  fracture  mechanics.  Crack  propagation  and 
fatigue  analyses  of  lap  joints  under  bending  and  tensile  loads  are  presented  in  the  third 
chapter  of  the  report. 

1.4  Contributions  of  Research 

The  specific  contributions  of  the  present  research  are: 

1 .  A  numerical  model  is  presented  to  compute  the  static  failure  strength  of  laser- 
welded  lap  joints.  The  stress  intensity  factors  at  the  sharp  notches  between  the 
welded  plates  are  determined  using  linear  elastic  fracture  mechanics  principles. 

2.  Parametric  studies  of  the  fracture  strength  of  several  different  laser-welded  lap 
configurations  are  conducted  and  the  results  are  used  to  develop  a  metamodel  for 
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the  rapid  analysis  of  lap  joints  based  on  its  geometric  dimensions.  A 
methodology  is  presented  for  designing  lap  joints  under  prescribed  loads  using  the 
metamodel. 

3.  Crack  propagation  and  fatigue  analysis  of  lap  joint  configurations  is  performed 
using  a  finite  element  program  called  FRANC2D.  Parametric  studies  of  the  crack 
trajectories  have  been  conducted  for  different  weld-widths  and  loading  conditions. 

4.  The  fatigue  life  of  laser-welded  lap  joints  has  been  analyzed  using  the  Paris  Law. 
It  is  demonstrated  that  the  weld-width  plays  a  significant  role  on  the  fatigue 
strength  of  the  lap  joint. 
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2.  ANALYSES  OF  LAP  JOINT  CONFIGURATIONS 


2.1  Introduction 

In  this  section  a  metamodel  is  presented  for  the  evaluation  of  the  fracture  strength 
of  a  laser  welded  lap  joint  of  arbitrary  dimensions  without  having  to  perform  a  dedicated 
finite  element  analysis  for  the  geometry  in  question  every  time.  This  is  achieved  by 
performing  detailed  finite  element  analyses  of  a  limited  number  of  joints  of  varying 
dimensions  and  utilizing  the  results  to  develop  a  higher  dimensional  curve  fit  using  radial 
basis  functions. 

It  is  difficult  to  obtain  a  closed-form  analytical  expression  for  the  fracture  strength 
of  the  laser  welded  lap  joints  considered  here.  A  finite  element  model  is  used  to  compute 
the  stress  intensity  factors  at  the  notches  and  determine  the  fracture  strength  of  the  joint. 
First,  convergence  studies  are  performed  on  representative  model  to  establish 
requirements  for  obtaining  accurate  results  from  the  finite  element  analyses.  The 
dependence  of  stress  intensity  factor  on  joint  dimensions  and  load-type  is  then 
investigated  by  performing  parametric  studies  using  finite  element  models.  To 
compensate  for  the  lack  of  an  analytical  expression,  a  metamodel  is  formed  to  interpolate 
the  results  obtained  from  the  parametric  studies.  The  metamodel  is  optimized  using  a 
single-objective  genetic  algorithm  and  validated  by  finite  element  results  evaluated  at 
sample  joint  dimensions  which  were  not  used  in  construction  of  the  metamodel. 

2.2  Problem  Formulation 

Fig.  2.1  depicts  the  type  of  connection  geometry  under  consideration  which  is 
referred  to  as  a  simple  lap  joint.  The  dimensional  parameters  here  are  welded  length  (w), 
non-welded  length  (a)  and  the  thicknesses  of  the  top  (h)  and  bottom  plates  (//).  The  width 
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of  the  joint  in  the  ^-direction  is  assumed  to  be  much  larger  than  the  plate  thickness  and 
the  loads  are  assumed  to  be  independent  of  the  ^-coordinate.  The  finite  element  analyses 
will  be  based  on  linear  elasticity.  Frictionless  hard  contact  is  assumed  over  the  non- 
welded  segments  of  the  plates.  It  is  anticipated  that  the  fracture  strength  will  have  a 
nonlinear  dependence  on  the  applied  load  due  to  the  contact  conditions.  The  material  of 
the  connection  plates  is  chosen  to  be  steel  as  previously  mentioned. 

2.3  Stress  Analysis  of  Laser  Welded  Lap  Joints 

A  detailed  analysis  of  a  representative  laser  welded  lap  joint  is  performed  using 
the  finite  element  program  ABAQUS,  to  quantify  the  stresses  in  the  welded  region.  The 
element  type  used  is  CPE4R,  a  four-noded  element  which  is  discussed  in  further  detail  in 
Section  2.6.2.  The  dimensions  of  the  representative  joint  are  chosen  as  h  =  H  =  w  =  a  = 
1".  Other  joint  configurations  will  be  considered  later  in  this  chapter.  The  plates  are 
chosen  to  be  sufficiently  long  in  the  *i -direction  in  order  to  minimize  the  influence  of  the 
edge  effects  on  the  stresses  in  the  welded  region.  The  entire  structure  is  subjected  to  an 
applied  counter-clockwise  moment  of  M  +  =  1  lb-in/in,  direction  of  which  is  shown  in 
Fig.  2.2(a).  Hard  contact  conditions  are  enforced  in  the  non-welded  regions  between  the 
two  welded  plates.  This  is  accomplished  by  specifying  contact  surfaces  on 
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Figure  2.1-  Laser  welded  lap  joint  configuration. 
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Figure  2.2:  Stress  analysis  of  a  laser  welded  lap  joint  under  an  applied  moment:  a)  loading  configuration,  b)  transverse 

normal  stress  and  c)  transverse  shear  stress. 


the  two  sides  of  the  interface,  and  thus,  penetration  of  one  surface  into  another  is 
prevented.  A  very  fine  mesh  is  used  to  ensure  that  the  results  have  converged.  Fig.  2.2(b) 
and  (c)  depict  the  transverse  normal  and  shear  stresses,  respectively,  at  the  interface 
between  the  plates  as  a  function  of  the  overlap  coordinates  X\.  It  is  observed  that  the 
stresses  are  singular  at  the  edges  of  the  weld  due  to  the  sharp  notches  created  by  the 
nature  of  the  geometry.  The  sharp  notches  induce  large  local  stresses  even  at  low  applied 
loads.  These  observations  justify  the  use  of  fracture  mechanics  concepts  for  this  joint 
analysis. 

2.4  Fracture  Criterion 

The  fracture  behavior  of  a  lap  joint  is  determined  by  the  stress  intensity  factors  at 
the  notches  located  at  the  left  and  right  edges  of  the  laser  welded  region.  The  high 
transverse  normal  stress  022  and  transverse  shear  stress  an  at  those  locations  may  result 
in  fracture  which  is  governed  by  both  mode-I  and  II  stress  intensity  factors  ( K )  and  Kn).  A 
combination  of  fracture  modes  is  referred  to  as  a  mixed  mode.  Under  mode-I  conditions, 
growth  occurs  in  the  direction  of  the  crack.  Under  mixed  mode  conditions,  experimental 
observations  indicate  that  both  Kj  and  Kn  play  a  major  role  in  crack  growth  and  that  the 
crack  does  not  grow  in  its  original  direction.  Both  energy  and  near  tip  stress  approaches 
are  commonly  used  in  the  study  of  mixed  mode  fracture.  In  the  energy  approach,  the 
crack  is  assumed  to  grow  in  the  direction  that  maximizes  the  energy  release  rate.  This 
approach  becomes  very  tedious  for  crack  branching  under  mixed  mode  fracture.  Several 
failure  theories  have  been  proposed  for  modeling  mixed  mode  fracture  that  provides 
results  consistent  with  experimental  observations.  These  include: 
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•  Elliptical  model. 

•  Maximum  tensile  stress  criterion  (MS-criterion). 

•  Strain  energy  density  criterion  (S-criterion). 

•  Maximum  energy  released  rate  criterion  (ME-criterion). 

We  utilize  the  maximum  tensile  stress  criteria  proposed  by  Erdogan  and  Sih  (1963) 
which  assumes  that  (i)  crack  extension  occurs  in  the  direction  at  which  the 
circumferential  stress  Gee  takes  the  maximum  value  with  respect  to  9  and  (ii)  fracture 
takes  place  when  the  intensity  of  the  maximum  circumferential  stress  Gee  is  equal  to  a 
value  that  leads  to  mode-I  fracture.  For  mixed  mode  fracture  the  stresses  near  the  crack 
tip  are 
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where  r  and  0  are  the  polar  coordinates  with  respect  to  the  crack  tip  of  the  material  point, 
as  shown  in  Fig.  1.1.  Using  the  direction  cosine  matrix,  the  stresses  can  be  transformed 
from  Cartesian  to  polar  coordinates  resulting  in 


G  ee  ~ 


°r0  = 


1  9 

- - cos  — 

Jim  2 

1  9 

, - cos— 

Jim  2 


2  9  3 


K,  cos - K„  sin  9 

2  2 


^KIsm0  +  ^KIl(3cos0-l) 


(4) 


From  equation  (4)  we  deduce  that 


doge 

dd 


(5) 
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This  shows  that  at  some  value  6  =  d0,  <Joe  reaches  its  maximum  when  -^7  =  are  -  0.  The 

do 


crack  initiation  angle  60  will  be  such  that  are  =  0,  and  thus,  from  the  second  equation  in 
(4)  the  following  condition  is  obtained  for  0o : 

K,  sin0o  +  KIt(3cos0o  - 1)  =  0.  (6) 

According  to  the  assumption  of  Erdogan  and  Sih  (1963),  fracture  takes  place  when 
(cr^)max  becomes  K/cN2jrr.  This  yields  the  fracture  criterion 


K.  cos2  —  -  —  K„  sin  0o  - - — - . 

2  2"  cos(0o  /  2) 


(7) 


Equation  (7)  represents  the  mixed  mode  fracture  criterion  where  K/c  is  the  critical 
mode-I  stress  intensity  factor  and  60  is  obtained  from  the  auxiliary  equation  (6).  It  is 
noted  that  only  one  material  property  K/c  is  involved  in  the  criterion.  If  the  first  of  the  two 
extremes  of  0o  is  computed,  this  being  associated  with  mode-I  or  when  Ku=  0,  then  0o=  0 
is  obtained  implying  that  the  crack  extends  in  its  original  direction.  For  the  other  extreme 
when  pure  mode-II  conditions  exist  (Kj  =  0),  the  computed  crack  extension  direction  is 
60  =  -70.5°.  Thus  it  can  be  inferred  that  for  a  mixed  mode  condition,  the  crack  initiation 
angle  will  be  between  -70.5°  and  0°. 

2.4.1  Critical  Load  Calculations 

No  fracture  mechanics  study  is  complete  without  discussing  the  critical  load  for 
fracture  initiation.  Under  linear  elastic  fracture  mechanics  theory,  it  can  be  said  that  the 
stress  intensity  factors  are  directly  proportional  to  the  applied  load.  This  assumes  that  no 
nonlinearity  is  present  in  the  structure.  The  stress  intensity  factors  K/Mo  and  KuMo  are 
obtained  via  analysis  of  the  structure  when  a  moment  M0  is  applied  at  the  free  end  of  the 
upper  plate  and  the  other  end  is  restrained.  According  to  linearity,  the  mode-I  stress 
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intensity  factor  for  a  moment  Mc  is  Mc-KiMo/M0.  By  substituting  this  expression  for  stress 
intensity  factor  in  equation  (7),  the  critical  moment  can  be  computed  as: 


Mc  = 


K 


lc 


cos- 


0 


k; 


Mn 


.cos 


0 '  3  K 


II 


2  2  M„ 


.sintf 


(8) 


This  critical  moment  represents  the  failure  load  of  the  structure.  This  expression  is  very 
useful  for  designing  joints  so  that  a  desired  factor  of  safety  for  a  fracture  type  failure  can 
be  achieved. 


2.5  Stress  Intensity  Factor  for  a  Finite  Specimen 

The  stress  intensity  factor  solution  for  some  basic  infinite  media  is  presented  in 
equation  (2).  The  solution  presented  there  is  valid  when  the  crack  length  is  much  smaller 
than  the  in-plane  dimensions  of  the  specimen.  If  the  crack  length  is  similar  in  magnitude 
to  characteristic  specimen  dimensions,  however,  the  stress  intensity  factor  becomes 
influenced  by  the  geometry  and  size  of  the  part.  In  these  cases  analytical  results  are 
difficult  to  obtain  and  thus  numerical  methods  (e.g.  the  finite  element  method)  are 
employed  to  calculate  stress  intensity  factors  at  the  cracks  in  finite  size  media. 


2.5.1  Stress  Intensity  Factor  by  the  Finite  Element  Method 

The  finite  element  method  has  become  a  useful  technique  for  solving  fracture 
mechanics  problems  of  various  types  including  those  of  the  linear  elastic  or  elasto-plastic 
variety.  Finite  element  approaches  to  fracture  mechanics  problems  can  be  described  as 
either  direct  or  indirect.  In  the  direct  method,  stress  intensity  factors  are  calculated 
directly  from  the  solution,  while  in  the  indirect  method,  the  strain  energy  release  rate  is 
calculated  and  stress  intensity  factors  are  inferred  from  it.  Among  the  two,  the  indirect 
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technique  is  the  most  popular  and  efficient.  There  are  several  methods  discussed  in  the 
literature  for  evaluating  the  energy  release  rate  for  finite  element  models.  These  methods 
are  based  on  changes  in  potential  energy,  compliance  or  stiffness  for  the  crack  before  and 
after  crack  extension.  Irwin  created  a  technique  called  the  Crack  Closure  Integral  in 
which  he  defined  the  energy  released  as 

da  i 

G,da  =  2\-ayyuydx.  (9) 

o  z 

The  above  expression  simply  states  that  the  work  done  by  Oyy  (the  far  field 
stresses  in  the  transverse  direction)  in  closing  the  extended  crack  length  da  by  the 
displacement  uy  is  equal  to  the  energy  released,  Gj.da.  A  graphical  depiction  of  these 
quantities  is  shown  in  Fig.  2.3.  The  factor  2  in  the  preceding  expression  signifies  that 
both  the  crack  surfaces  are  taken  into  account. 

This  method  is  not  widely  used,  perhaps  because  it  required  two  separate  analyses 
to  determine  the  stress  intensity  factor  K/.  Regardless,  this  method  is  based  on  the 
assumption  that  the  energy  A E  released  when  the  crack  is  extended  from  a  to  a+Aa,  is 
equal  to  the  energy  required  to  close  the  crack  between  location  /  and  i  shown  in  Fig. 
2.4(a).  For  a  crack  modeled  with  two-dimensional,  four-noded  elements  as  shown  in  Fig. 
2.4(a),  the  work  AE  required  to  close  the  crack  along  one  element  side  can  be  expressed 
as 

AE  =  ±[FyAvl  +  FxAut\  (10) 

where  Fx  and  Fy  are  the  shear  and  opening  forces  at  the  crack-tip  node  /  and  Av1  and  A  ul 
are  the  differences  in  the  shear  and  opening  displacements  at  node  /  shown  in  Fig.  2.4(b). 
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Figure  2.3:  Depiction  of  crack  closure  integral. 
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Figure  2.4:  Irwin’s  Crack  Closure  Technique,  a)  step-I  and  b)  step-II. 
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In  the  first  analysis,  a  finite  element  model  of  the  cracked  structure  is  created  with 
the  desired  loads  and  boundary  conditions,  and  the  forces  Fx  and  Fy  are  computed  at  the 
crack-tip  /,  shown  in  Fig.  2.4  (a).  In  the  second  analysis,  the  model  is  re-run  with  the 
crack  being  extended  to  a+Aa  and  displacement  of  point  /  is  noted.  It  follows  that  if  the 
forces  Fx  and  Fy  were  not  present,  then  point  /  would  displace  by  the  amounts  Av7  and  A  u‘ 
in  the  vertical  and  horizontal  directions  respectively. 

To  eliminate  the  need  for  multiple  analyses  the  crack  closure  technique  is 
modified  in  such  a  fashion  that  only  one  analysis  is  required  to  obtain  values  of  both  K/ 
and  K/i.  This  method  is  called  the  Modified  Crack  Closure  Technique  and  it  also  falls  into 
the  category  of  indirect  methods.  The  method  does  not  require  the  calculation  of  stresses 
and  conventional  constant  strain  elements  can  be  used  successfully. 

In  Fig.  2.4(a),  if  Fy  and  Fj  are  vertical  and  horizontal  forces  at  node  /  before 
extension  and  u1,  u' ,  v7  and  v7  are  the  respective  horizontal  and  vertical  displacements  of 
nodes  /  and  /’  after  crack  extension  (Fig.  2.4(b)),  then  equation  (10)  can  be  manipulated 
to  express  the  strain  energy  release  rate  in  terms  of  nodal  quantities,  giving 


1 


2A a  y 
1 


2A a 


F'x  (u1  -  u1 ). 


(11) 


The  superscripts  in  this  equation  denote  the  node  location  as  shown  in  Fig.  2.4.  For  the 
modified  crack  closure  technique,  a  very  small  Aa  is  required  and  should  be  on  the  order 
of  Aa/a  <  0.05.  Rybicki  and  Kanninen  (1977)  assumed  that  if  Aa/a  is  small  and  Aa  kept 
the  same  before  and  after  the  crack-tip,  the  nodal  displacements  at  /  and  /’  after  crack 
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extension  can  be  approximated  by  the  nodal  displacements  at  j  and  j’  before  the  crack 
extension.  Thus,  equation  (11)  can  now  be  written  as 


(12) 


With  this  method,  only  a  single  analysis  needs  to  be  run  in  order  to  obtain  K/  and  A"//  for 


any  crack  in  the  structure. 


2.5.1. 1  Correction  for  Unequal  Element  Lengths  at  the  Crack  Tip 

While  evaluating  Fyl  and  Fj  it  is  convenient  to  maintain  //  =  Aa  and  also  to  keep 


elements  (1),  (2),  (3)  and  (4)  the  same  size  (Fig.  2.5).  It  is  not  difficult  to  maintain  equal 
length  elements  around  the  crack  tip  by  controlling  the  finite  element  mesh.  One  might 
do  this  by  simply  creating  a  square  mesh  partition  around  the  crack  tip.  However,  when 
unequal  length  elements  cannot  be  avoided,  a  modification  to  equation  (12)  allows  one  to 
use  elements  of  differing  lengths  around  the  crack  tip.  Referring  to  Fig.  2.5,  F 'x  and  Fv  are 
the  required  nodal  forces  at  i  for  //  =  Aa,  but  for  the  case  of  unequal  element  lengths  this 
is  not  exactly  true  and  Fx  and  Fy  are  simply  the  nodal  forces  obtained  from  the  finite 
element  model.  To  obtain  the  correct  required  nodal  forces,  the  following  correction 
factors  can  be  used  (Rybicki  and  Kanninen,  1977): 


(13) 
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Figure  2.5:  Correction  method  in  the  Modified  Crack  Closure  Technique  for  unequal 

length  elements. 
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Substituting  these  expressions  into  equation  (12)  gives  the  required  result,  which  is 


(14) 


2.5.1. 2  Relationship  between  G  and  K  by  Crack  Closure  Method 

As  stated  by  Irwin,  if  a  crack  extends  by  a  small  amount  Aa,  the  change  in  energy 
in  the  process  is  equal  to  the  work  required  to  close  the  crack  to  its  original  length.  This 
can  be  expressed  as 


(15) 


where  ay  and  rxy  are  the  tensile  and  shear  stresses,  respectively,  at  node  /  before  extension 
and  u  and  v  are  sliding  and  opening  displacements  at  node  /  after  the  crack  extension 
(Fig.  2.4(a)).  The  first  and  second  parts  of  the  total  energy  release  rate  G  are  essentially 
two  separate  contributions  for  mode-I  and  mode-II,  respectively.  Noting  that 
Westergaard’s  expressions  for  oy,  txy,  u  and  v  are 


and  then  substituting  these  expressions  into  equation  (15)  will  give 


(16) 
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rc  =  < 


(16-a) 


3  -  4v  plane  strain 
3  -v 

-  plane  stress 

[\  +  v 

2.5.2  Stress  Intensity  Factor  and  its  Dependence  on  Geometry 

Obtaining  closed  form  solutions  for  the  stress  intensity  factors  K/  and  Ku  for  the 
lap  joint  configuration  would  entail  the  solution  of  complex  two-dimensional  elasticity 
problem.  It  is  particularly  difficult  to  obtain  an  analytical  solution  to  the  governing  partial 
differential  equations  due  to  the  hard  contact  conditions  at  the  interface  of  the  non- 
welded  regions  of  the  joint.  In  the  present  work,  the  finite  element  package  ABAQUS  (v. 
6.6)  is  utilized  to  evaluate  the  stresses  and  the  stress  intensity  factors  at  the  crack  tips.  A 
two-dimensional  model  with  plane  strain  conditions  is  analyzed  with  a  very  fine  mesh 
consisting  of  approximately  300,000  elements.  The  non-welded  interface  between  the 
plates  is  modeled  with  hard  normal  contact  conditions  with  negligible  friction.  The  left 
end  of  the  lower  plate  is  held  fixed. 

The  stress  intensity  factors  at  the  notch  tips  are  computed  using  the  Modified 
Crack  Closure  Technique  and  the  finite  element  method.  In  the  MCC  technique,  a 
regular  mesh  of  four-noded  quadrilateral  elements  is  used  around  the  crack  tip. 
According  to  equations  (12)  and  (16),  the  stress  intensity  factors  are  computed  from  the 
nodal  forces  at  the  notch  tip  and  the  relative  displacement  of  the  nodes  adjacent  to  the 
notch  tip. 

By  executing  several  analyses  with  different  joint  dimensions  and  load  types  it  is 
found  that  the  stress  intensity  factors  at  the  notch  tips  vary  with  the  change  in  the 
dimensional  parameters  (w,  a,  H  and  h)  of  the  joint.  This  fact  gives  credence  to  the 
statement  in  Section  2.4  that  stress  intensity  factors  are  functions  of  the  geometric 
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dimensions  of  the  lap  joint  and  the  nature  of  the  applied  loads.  To  aid  in  determining  the 
dependence  of  the  stress  intensity  factor  on  joint  dimensions,  the  four  aforementioned 
parameters  are  normalized  with  respect  to  h,  thus  reducing  the  number  of  variables  from 
four  to  three.  Similar  to  /?  in  equation  (2),  a  geometric  configuration  factor  /  is  defined 
which  depends  on  the  three  non-dimensional  geometric  variables, 

r  M+ ,1  /  //  _>v  a\ 

J i  V  /»’/»»  A/*  (17) 


where  the  subscript  i  =  /,  II  denotes  the  fracture  mode.  The  first  superscript  of  the 
function  f  indicates  the  direction  of  the  applied  moment  which  can  either  be  M +  for  a 
positive  moment  (counter-clockwise)  or  M  '  for  a  negative  moment  (clockwise).  The 
second  superscript  of  the  function  f  indicates  the  location  of  the  notch  which  can  either 
be  /  for  the  notch  on  the  left  edge  of  the  weld  or  r  for  the  notch  on  the  right  edge.  In  a 
work  by  Suo  and  Hutchinson  (1990),  an  analytical  solution  for  the  stress  intensity  factors 
for  a  semi-infinite  crack  between  two  infinite  elastic  layers  was  obtained  that  provided 
some  insight  as  to  how  to  properly  construct  a  mathematical  model  for  defining  K/  and 
K//.  Based  on  the  functional  form  of  the  stress  intensity  factors  obtained  by  Suo  and 
Hutchinson  (1990)  and  other  dimensional  considerations,  the  following  expressions  for 
the  stress  intensity  factors  are  assumed: 


M  M+,i  H  w  q 

3/2  Ji  V  ,  5  i  5  , 


h  h  h 


(18) 


where  M  is  the  magnitude  of  the  applied  moment  and  f  are  previously  defined  geometric 
configuration  factors.  In  the  present  work,  the  functions  f  are  evaluated  using  the  finite 
element  method  at  discrete  values  of  the  geometric  parameters.  The  values  obtained  from 
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the  parametric  studies  are  later  used  to  construct  a  three-dimensional  metamodel  to 
predict  stress  intensity  factors  for  arbitrary  joint  dimensions. 

2.5.3.  Critical  Load  Scenario  for  a  Fixed  Joint  Dimension 

The  objective  of  this  study  is  to  check  if  LEFM  is  valid  for  our  mathematical 
model.  To  do  that,  one  needs  to  know  the  plasticity  effects  in  a  specific  joint  due  to 
application  of  the  critical  moment.  In  order  to  obtain  a  converged  critical  moment  of  a 
joint,  a  convergence  study  is  performed.  A  particular  joint  geometry  (H/h  =  1,  w/h  =  0.6, 
a/h  =  1)  is  selected  with  upper  plate  thickness  h  =  1”,  and  critical  moment  is  calculated 
utilizing  the  stress  intensity  factors  obtained  by  applying  the  unit  positive  moment.  The 
finite  element  model  with  boundary  conditions  and  load  is  shown  in  Fig.  2.6.  The  critical 
moment  thus  obtained  is  applied  to  the  joint  and  a  new  critical  moment  is  evaluated  based 
on  new  stress  intensities.  The  process  is  iterated  until  convergence  is  established  and 
variation  in  the  results  reaches  less  than  0.01%.  Data  corresponding  to  this  study  is 
tabulated  and  shown  in  Table  2.1.  A  graphical  depiction  of  the  convergence  results  is 
shown  in  Fig.  2.7.  According  to  the  assumed  linear  relation  between  stress  intensity 
factor  and  applied  load  (shown  in  equation  (18)),  if  the  stress  intensity  factors  obtained  at 
unit  load  be  extrapolated,  one  should  get  very  close  values  to  the  stress  intensity  factors 
obtained  at  critical  load.  However,  these  values  would  not  be  close  due  to  the  non¬ 
linearity  enforced  due  to  contacts. 
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Figure  2.5:  Convergence  of  the  critical  moment  for  the  joint. 


is  on  finding  the  amount  of  load  which  should  be  applied  to  the  joint  which  gives  lowest 
extrapolated  error  but  would  still  follow  LEFM  principles.  If  the  stress  intensities  at  an 
initial  moment  of  M0  =  1  lb-in/in  are  extrapolated,  it  is  seen  that  the  actual  stress 
intensities  for  mode-I  at  the  true  critical  load  are  4.9%  higher  than  those  obtained  from 
the  extrapolation.  This  explains  that  for  this  particular  geometry,  extrapolating  the  results 
of  a  unit  applied  moment  gives  underestimated  fracture  calculations.  However,  if  the 
analysis  is  re-run  with  higher  load,  the  extrapolated  stress  intensity  factors  come  closer  to 
the  actual  values  obtained  at  that  load.  A  set  of  such  analyses  have  been  performed  on  the 
same  geometry  with  increasing  moments  up  to  the  critical  moment  of  the  joint,  to 
recognize  the  appropriate  moment  to  be  applied  for  better  predictions.  The  results  of  this 
study  are  shown  in  Tables  2.2(a)  and  (b)  for  Kj  and  Kn  respectively. 

As  seen  in  the  tables,  each  load  step  is  selected  one  by  one  and  stress  intensities 
are  extrapolated  based  on  the  selected  load.  Errors  from  the  true  stress  intensity  values 
are  computed  in  each  increment  case.  With  these  results  in  hand,  it  is  then  at  the 
discretion  of  the  analyst  as  to  how  much  error  can  be  tolerated.  Selecting  the  moment  of 
20,000  lb-in/in,  gives  1 .9%  underestimated  and  0.7%  overestimated  values  for  K/  and  Kn 
respectively.  Prior  to  moving  forward,  one  must  check  to  see  if  LEFM  is  still  applicable. 
Von  Mises  stress  contour  plots  have  been  obtained  for  4147  steel  in  order  to  check  the 
plastic  zone  size  at  around  55%  of  the  critical  moment  (20,000  lb-in/in)  of  the  joint. 
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the  plastic  zone  is  larger  (0.059”)  than  the  left  lobe  (0.0243”)  due  to  the  apparent  stress 
distribution.  It  is  quite  obvious  but  worth  mentioning  that  the  plastic  zone  size  will 
increase  for  a  more  ductile  material. 

Gdoutos  (2005)  has  derived  a  relation  for  the  plastic  zone  radius  as  a  function  of  8 
which  is  in  the  form 


2ncr, 


K)  cos2 1|(1  -  2v)2  +  3  sin2 1|  +  K,KU  sin  8^,  cos  8  -  (1  -  2vf}+ 


2  @  //i  ->  \2  «  y  O' 


j  3  +  sin  — u  1  ~  2v)  -9cos  — 


.(19) 


This  relation  is  valid  for  mixed  mode,  plane  strain  conditions.  Fig.  2.9  depicts  the 
shape  of  the  theoretical  plastic  zone  in  polar  coordinates  for  the  stress  intensities  obtained 
at  an  applied  moment  of  20,000  lb-in/in.  As  expected,  the  zone  size  changes  for  different 
material  strengths  as  evinced  by  the  fact  that  rp  depends  on  the  yield  strength  (<jy).  The 
theoretical  plastic  zone  size  and  shape  is  compared  with  the  one  shown  in  the  contour 
plots  from  the  joint  analysis.  The  gross  dimensions  of  the  plastic  zone  computed  from  the 
analytical  and  finite  element  models  are  relatively  similar. 
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Figure  2.8:  Contour  plot  showing  gray  plastic  zone  for  4147  steel  of  yield  strength  137ksi. 
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Figure  2.9:  Theoretical  plastic  zone  shape  and  size  for:  a)4147  steel 

and  b)HSLA-65  steel. 
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Since  the  plastic  zone  size  is  seen  within  l/50th  of  the  thickness  dimensions  of  the  joint, 
one  can  safely  use  LEFM  for  4147  steel. 

2.6  Finite  Element  Modeling 
2.6.1  Analysis 

A  commercially  available  finite  element  package  ABAQUS  (2005),  is  used  to 
perform  numerical  modeling.  Some  of  the  major  modeling  details  are  provided  in  this 
section  to  provide  information  so  that  this  work  can  be  duplicated  using  other  finite 
element  packages.  A  two-dimensional  model  is  prepared  by  constructing  two  separate 
plates  and  joining  them  in  the  middle  of  the  overlap  to  achieve  the  desired  joint 
configuration,  as  already  shown  in  Fig.  2.1.  To  simulate  the  weld,  the  nodes  on  the 
welded  length  of  the  overlapped  plates  are  assigned  the  same  degrees  of  freedom  (Fig. 
2.10).  The  non-welded  overlapped  regions  have  been  defined  as  frictionless  hard  surfaces 
to  ensure  non-penetrating  contact.  Both  the  top  and  the  bottom  plates  are  assumed  to  be 
of  same  isotropic  material,  steel,  with  Young’s  modulus  and  Poisson’s  ratio  of  29.7xl03 
ksi  and  0.29  respectively.  The  finite  element  model  is  fixed  at  the  left  edge  of  the  lower 
plate  and  point  loads  are  applied  at  the  right  free  edge  of  the  upper  plate  to  simulate  a 
couple  (refer  Fig.  2.6) 

A  combination  of  triangular  (continuum,  plane-strain,  3-noded-CPE3,  away  from 
the  notches)  and  quadrilateral  (continuum,  plane-strain,  4-noded,  reduced  integration, 
hourglass  control-CPE4R,  near  the  notch  tips)  elements  are  used,  each  possessing  linear 
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Figure  2.10:  Deformed  finite  element  model  of  the  joint  (H/h=l,  w/h=0.6  and  a/h=l). 
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shape  functions.  This  combination  is  used  to  reduce  the  number  of  elements  in  each 
model  and  reduce  the  analysis  runtime.  The  selection  of  elements,  their  number  and  their 
sizes  will  be  discussed  in  detail  in  Section  2.6.2  and  2.6.3.  The  geometry  and  properties 
are  defined  and  the  model  is  set  to  run  under  plane  strain  conditions.  Concerning  the 
solution  method,  ABAQUS/Standard  provides  its  own  algorithm  of  automatic 
incrementation  instead  of  asking  the  user  to  select  among  Newton  or  Modified  Newton’s 
methods.  This  default  setting  provides  reliable  results  and  a  fast  convergence  rate.  A 
direct  matrix  solver  is  used  at  each  iteration  which  reduces  the  solution  time  and  space 
requirements  by  taking  maximum  advantage  of  matrix  sparsity.  It  is  important  that  the 
appropriate  output  options  are  selected  so  that  nodal  forces  and  displacements  can  be 
retrieved  at  the  crack  tips. 

2.6.2  Element  Type 

ABAQUS  provides  several  different  types  of  two-dimensional  elements.  For 
structural  applications  these  include  plane  stress  elements  and  plane  strain  elements. 
Plane  stress  elements  can  be  used  when  the  thickness  of  a  body  or  domain  is  small 
relative  to  its  lateral  (in-plane)  dimensions.  Plane  strain  elements  can  be  used  when  it  is 
safe  to  assume  that  the  strains  in  a  loaded  body  or  domain  are  functions  of  planar 
coordinates  only,  and  that  the  out-of-plane  normal  and  shear  strains  are  equal  to  zero. 
Plane  strain  elements  must  be  defined  in  the  X-Y  plane  and  all  loading  and  deformation 
is  restricted  to  this  plane  as  well.  This  modeling  method  is  generally  used  for  bodies  that 
are  very  thick  relative  to  their  lateral  dimensions,  like  the  model  under  consideration. 

There  are  five  factors  to  characterize  the  behavior  of  an  element.  These  are 
namely  the  element  family,  the  number  of  degrees  of  freedom,  the  number  of  nodes,  the 
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element  formulation  type  and  the  integration  procedure  used.  Fig.  2.1 1  shows  the  element 
families  available  in  ABAQUS  that  are  most  commonly  used  in  a  stress  analysis.  For  our 
analysis,  where  the  plates  are  considered  to  be  of  infinite  width  and  it  is  important  to 
model  contacts,  two-dimensional  continuum  elements  CPE4R  and  CPE3  are  selected, 
where  the  first  letter  or  letters  indicate  the  element  family.  For  example,  CPE4R  is  a 
continuum,  plane  strain  element  (C  -  continuum,  PE  -  plane-strain).  The  degrees  of 
freedom  are  the  fundamental  variables  calculated  during  the  analysis.  For  a 
displacement/stress  element  the  degrees  of  freedom  are  the  translations.  At  any  point  in 
the  element  other  than  the  nodes,  the  displacements  or  other  degrees  of  freedom  are 
obtained  by  interpolating  from  the  nodal  degrees  of  freedom.  Usually  the  interpolation 
order  is  determined  by  the  number  of  nodes  used  in  the  element.  Elements  that  have 
nodes  only  at  their  comers,  such  as  a  rectangular  4-node  2D  continuum  element,  use 
linear  interpolation  in  each  direction  and  are  labeled  linear  elements  or  first-order 
elements. 
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The  formulation  of  an  element  refers  to  the  mathematical  theory  used  to  define 
the  element's  behavior.  Since  adaptive  meshing  techniques  are  not  used  in  our  work,  all 
of  the  displacement/stress  elements  used  here  are  based  on  a  Lagrangian  formulation,  i.e. 
the  element  deforms  with  the  material.  With  this  type  of  formulation  large  scale 
deformations  cannot  be  modeled  correctly. 

The  finite  element  program  ABAQUS  uses  numerical  techniques  to  integrate 
various  quantities  over  the  volume  of  each  element,  thus  allowing  complete  generality  in 
material  behavior  (Section  3.2.1,  ABAQUS  Theory  Manual).  Using  Gaussian  quadrature 
for  most  elements,  including  CPE4R  and  CPE3,  ABAQUS  evaluates  the  material 
response  at  each  integration  point  in  each  element.  Some  continuum  elements,  such  as 
CPE4R,  have  the  ability  to  use  either  full  or  reduced  integration.  Reduced  integration 
decreases  the  number  of  constraints  introduced  by  an  element  when  there  are  internal 
constraints  in  the  continuum  theory  being  modeled,  such  as  incompressibility  or  the 
Kirchhoff  transverse  shear  constraints  in  the  case  of  solid  elements  used  for  analyzing 
bending  problems.  In  such  applications  fully  integrated  elements  will  “lock”  and  will 
exhibit  a  very  stiff  response,  thus  giving  very  poor  results.  In  these  cases,  the  reduced 
integration  version  of  an  element  will  often  work  well.  The  letter  R  at  the  end  of  an 
ABAQUS  element  name  signifies  a  reduced  integration  element.  Although  reduced 
integration  lowers  the  cost  of  forming  an  element,  the  drawback  is  that  except  in  one¬ 
dimensional  and  in  axisymmetric  geometries  modeled  with  higher  than  first-order 
elements,  the  element  stiffness  matrix  will  be  rank  deficient.  This  often  creates  singular 
or  hourglass  modes  in  the  response.  When  reduced  integration  is  used  in  first-order 
elements  (4-node  quadrilateral),  hourglassing  often  makes  the  elements  unusable  unless 
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the  phenomena  is  controlled.  In  ABAQUS,  artificial  stiffness  and  damping  methods  are 
used  to  control  the  hourglass  modes  in  these  elements. 

The  CPE4R  and  CPE3  elements  are  applicable  for  many  applications.  These 
displacement/stress  elements  can  be  used  in  linear  or  complex  nonlinear  mechanical 
analyses  that  involve  contact,  plasticity,  and/or  large  deformations. 

2.6.3  Convergence  Study 

A  convergence  study  is  required  at  several  steps  of  the  process  to  verity  that  the 
mesh  being  used  is  sufficiently  refined.  The  first  study  is  performed  to  assess  the 
accuracy  of  the  Modified  Crack  Closure  Technique.  This  is  performed  using  two  separate 
center  crack  plate  models  subjected  to  the  remotely  applied  uniform  normal  and  shear 
tractions  for  K/  and  Kn  respectively.  The  dimensions  of  the  plate  in  the  model  are  made  to 
be  much  larger  than  the  crack  width  2a,  thus  allowing  the  use  of  analytic  solutions  for 
infinite  media  for  comparison  purposes.  The  element  types  and  analysis  procedure  used 
are  the  same  as  was  discussed  in  Section  2.6.1  and  2.6.2.  The  finite  element  solution  has 
been  compared  with  the  analytic  values  and  the  results  of  the  study  are  presented  in 
Tables  2.3(a)  and  (b).  The  errors  in  the  finite  element  results  are  plotted  as  a  function  of 
element  size  in  Fig.  2.12. 

Referring  to  Table  2.3(a),  a  mesh  size  of  Aa/a  =  0.0025  appears  to  give  good 
results  for  Kt,  as  the  reported  value  is  only  0.4%  off  from  the  analytic  value.  Further 
mesh  refinement  provides  only  diminishing  gains  in  accuracy.  Rather  than  decreasing  the 
element  size,  one  can  double  the  partition  area  close  to  the  notch  tip  while  keeping  the 
same  element  size  ( Aa/a  =  0.0025)  which  reduces  the  error  to  only  0.04%.  This 
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Error  (%) 


a)  Center  crack  convergence  study  for  Ki 


Ele.  Size(Aa) 

K, 

Analytical 

Error  % 

Aa/a 

0.003 

79.013 

79.461 

0.564 

0.015 

0.0009 

79.491 

79.461 

0.038 

0.0045 

0.0007 

80.09 

79.461 

0.792 

0.0035 

0.0005 

79.803 

79.461 

0.430 

0.0025 

b)  Center  crack  convergence  study  for  K» 


Ele.  Size(Aa) 

K„ 

Analytical 

Error  % 

Aa/a 

0.003 

78.446 

79.267 

1.036 

0.015 

0.0009 

78.912 

79.267 

0.448 

0.0045 

0.0007 

79.499 

79.267 

0.293 

0.0035 

0.0005 

79.221 

79.267 

0.058 

0.0025 

Table  2.3:  Convergence  study  for  the  center  crack  model  for:  a)  mode-1  stress  intensity 

factor  and  b)  mode-II  stress  intensity  factor. 


Aa/a 


Figure  2.12:  Convergence  study  of  stress  intensity  factors  for  the  center  crack  model. 
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essentially  means  defining  a  larger  region  around  the  crack  tip  for  which  there  is  a  finer 
mesh.  Table  2.3(b)  shows  the  convergence  for  the  mode-II  fracture  parameter,  and  for  an 
element  size  ratio  of  A  a/ a  =  0.0025  around  the  notch  tip,  the  error  is  only  0.06%. 

To  establish  if  the  small  element  size  and  mesh  partitions  used  in  the  comparison 
with  analytic  solutions  would  work  with  the  connection  geometry  under  consideration,  a 
similar  convergence  study  is  also  performed  on  a  selected  joint  configuration  (H/h  =  1 , 
w/h  =  1,  a/h  =  1).  Performing  a  convergence  study  on  the  model  with  element  size  ratios 
of  0.015,  0.0045,  0.0035,  0.0025  close  to  the  notch  tip,  yields  results  less  than  1%  off 
from  each  other.  Due  to  some  technical  issues  with  the  software’s  meshing  features,  it 
was  hard  to  obtain  an  element  size  ratio  around  the  crack  tip  any  smaller  than  0.0025. 
Normalizing  the  selected  element  size  ratio  in  terms  of  h,  Aa/h  =  0.000625  is  obtained.  In 
order  to  obtain  accurate  stress  intensity  factors  for  laser  welded  lap  joints,  an  even  finer 
mesh  is  used  in  subsequent  analyses  by  choosing  the  normalized  element  size  to  be 
Aa/h= 0.0005. 

2.6.4  Contact  Modeling 

In  general,  structural  surfaces  do  not  penetrate  each  other.  Thus,  in  order  to  obtain 
more  realistic  and  accurate  values  for  the  stress  intensity  factors,  it  is  important  to 
incorporate  contact  conditions  between  the  non-welded  regions  of  the  joint  in  the  finite 
element  model.  It  is  found  that  neglecting  contacts  may  provide  mode-I  stress  intensity 
factors  that  are  off  by  500-600%  from  the  true  values.  Typically,  the  mode-I  stress 
intensity  factor  value  is  lower  in  the  contact  model  since  the  normal  reaction  stresses  are 
distributed  over  a  wider  contact  region  in  the  joint.  When  contact  is  neglected,  all  stresses 
are  distributed  over  the  narrow  welded  region  which  results  in  large  stress  intensity 
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factors  at  the  notch  tips.  In  general,  the  mode-II  stress  intensity  factor  value  is  also 
smaller  in  the  contact  model,  though  the  difference  is  not  as  significant  as  it  is  for  the 
mode-I  stress  intensity  factor.  The  inclusion  of  contact  conditions  in  the  finite  element 
model  increases  the  computational  cost  due  to  the  iterative  solution  procedure  required  to 
solve  the  now  geometrically  nonlinear  problem.  In  addition,  the  inclusion  of  contact 
causes  a  somewhat  nonlinear  relationship  between  the  computed  stress  intensities  and  the 
load  applied.  This  discrepancy  was  discussed  in  detail  in  Section  2.5.3. 

If  a  linear  structure  is  subjected  to  different  types  of  loads,  then  it  is  possible  to 
evaluate  the  stress  intensity  induced  by  each  type  of  load  separately  and  then  sum  the 
solutions  to  get  the  total  stress  intensity  at  the  crack  tip.  However,  the  principle  of 
superposition  is  limited  to  only  linear  structures.  Even  though  our  structure  is  nonlinear,  a 
simple  test  is  performed  to  check  if  the  principle  of  superposition  can  be  applied  on  the 
joint  geometry  under  consideration,  thus  allowing  the  stress  intensity  factor  to  be 
calculated  for  various  loads  by  using  results  from  a  select  few  simple  cases. 

2.6.4.1  Superposition  of  Stress  Intensity  Factors 

Various  uniform  loading  cases  are  studied  in  literature,  but  in  actual  engineering 
applications,  a  body  is  often  subjected  to  non-uniformly  distributed  loads.  Stress  intensity 
factors  at  a  crack  in  the  body  under  arbitrary  loading  conditions  can  be  obtained  using  the 
method  of  superposition  for  linear  structures.  This  allows  the  evaluation  of  the  stress 
intensity  factor  of  a  desired  load  case  using  available  fundamental  solutions  for  bodies 
subjected  to  uniform  forces. 

In  an  attempt  to  model  the  real  life  loading  situations  faced  by  the  lap  joint 
studied  herein,  a  small  finite  element  modeling  experiment  is  performed.  It  is  assumed 
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that  to  simulate  the  response  of  the  joint  to  the  desired  variety  of  arbitrary  loads,  the 
results  from  only  two  types  of  fundamental  loads,  moment  and  axial  loads,  is  required. 
By  superposition  of  the  stress  intensity  factors  obtained  due  to  these  loads  separately,  one 
can  compute  the  stress  intensity  factors  at  the  notches  due  to  any  linear  combination  of 
these  loads.  Consider  a  lap  joint  configuration  fixed  at  the  left  end  of  the  lower  plate,  as 
shown  in  the  Fig.  2.13(a).  The  stress  intensities  in  this  joint  can  be  obtained  by 
superimposing  the  solutions  of  the  following  two  problems,  as  illustrated  in  Figs.  2.13(b) 
and  (c). 

A  joint  with  dimensions  of  H/h  =  2.5,  w/h  =  0.6,  and  a/h  =  1  is  selected.  Finite 
element  models  are  made  and  stress  intensity  factors  for  mode-I  and  II  are  computed  for 
the  left  and  right  notches,  for  both  the  loading  conditions  (a)  and  (b)  assuming  unit  tensile 
force  P,  with  contact  conditions  in  effect.  Table  2.4(a)  shows  the  results  obtained  for  K/ 
and  Kfi  values  at  the  left  and  right  notches  of  the  joint.  It  shows  that  there  is  a  large 
discrepancy  between  the  superimposed  and  directly  computed  mode-I  stress  intensities, 
the  reason  for  which  can  be  attributed  to  the  nonlinearity  created  by  the  contact 
conditions.  To  ascertain  this  reasoning,  the  contact  conditions  were  removed  from  all  the 
models  and  the  analyses  were  rerun.  Table  2.4(b)  displays  the  results  obtained  for  the 
linear  models.  The  results  verily  that  superposition  works  for  the  models  without 
contact.  This  study  proves  that  including  contact  phenomena  induces  geometric 
nonlinearity  in  the  connection.  Therefore,  the  principle  of  superposition  cannot  be 
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Figure  2.13:  Superposition  of  models,  for  computing  stress  intensity  factors  of  the  joint 

under  combined  loading. 
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(c) 

0.66 
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Left  end 

Right  end 

Contacts 

K, 

Kn 

K, 
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(b) 

0.112 

0.837 

0.134 

1.363 

(c) 

3.357 

0.325 

3.541 

-0.469 

Total 

3.469 

1.162 

3.675 

0.894 

Ja) _ 

3.459 

1.162 

3.404 

0.896 

Table  2.4:  Effect  of  contact  in  stress  intensity  factor  calculations:  a)  with  contact  b) 

without  contact. 
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utilized  in  this  work.  Unfortunately,  for  other  arbitrary  loading  conditions  of  interest,  one 
has  to  make  a  separate  model  for  computing  the  fracture  strength  of  the  connection. 


2.6.5  Parametric  Study 

Detailed  finite  element  analyses  are  performed  for  all  combinations  of  the  joint 
parameters  Hlh  =  [1,  2.5,  5],  w/h  =  [0.2,  0.6,  1]  and  alh  =  [0.25,  1,  2]  and  the  stress 
intensity  factors  are  evaluated  for  each  of  the  27  different  joint  configurations.  Table  2.5 
lists  the  normalized  mode-I  and  mode-II  stress  intensity  factors  at  the  left  and  right  notch 
for  all  three  ratios  of  Hlh  corresponding  to  a  positive  applied  moment.  Similar  results 
corresponding  to  a  negative  moment  are  shown  in  Table  2.6.  Tables  2.5  and  2.6  also  list 
the  mode  mixity  angle,  defined  as  tan‘'(Kn/K|),  for  various  joint  configurations.  Since 
the  critical  moment  is  a  dependent  parameter,  it  is  normalized  as 


Mc 


M  = 

c  h3l2.  K,c 

When  a  positive  or  negative  moment  is  applied  at  the  free  edge  of  the  upper  plate, 
one  can  visualize  one  side  of  the  overlap  opening  and  the  other  closing.  Cracks  that  close 
completely  gives  K/=  0  and  are  highlighted  as  gray  boxes  in  the  Tables  2.5  and  2.6.  The 
negative  stress  intensity  values  in  both  the  tables  for  a  closing  notch  signify  the  reverse 
stresses  of  shear.  For  the  M~  cases,  there  are  fewer  crack  closings.  By  zooming  in  on  the 
deformed  model  shown  in  Fig.  2.14,  this  behavior  can  be  explained  by  the  hinging  effect 
of  the  overhang  plate.  The  stresses  which  are  usually  carried  by  the  welded  length 
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101200 

1.322 

1 

0.6 

0.25 

4.53 

1.32 

21620 

0.282 

0.00 

-1.58 

59600 

1 

0.6 

1 

2.11 

1.30 

37080 

0.551 

0.00 

-1.30 

72600 

1 

0.6 

2 

1.66 

1.40 

40560 

0.701 

0.00 

-1.40 

67530 

1 

1 

0.25 

2.71 

1.40 

31040 

0.477 

0.00 

-1.59 

59560 

1 

1 

1 

1.85 

1.45 

37780 

0.666 

0.07 

-1.45 

63540 

1.523 

1 

1 

2 

1.63 

1.45 

40040 

0.728 

0.10 

-1.45 

62940 

1.499 

2.5 

0.2  0.25 

11.87 

1.21 

9047 

0.101 

0.00 

-1.43 

66200 

2.5 

0.2 

1 

3.28 

0.71 

31170 

0.212 

0.65 

-0.67 

91770 

0.801 

2.5 

0.2 

2 

2.43 

0.75 

39920 

0.299 

0.37 

-0.72 

103600 

1.100 

2.5 

0.6  0.25 

4.65 

1.18 

21579 

0.248 

0.02 

-1.12 

83343 

1.553 

2.5 

0.6 

1 

2.57 

1.15 

34380 

0.421 

0.00 

-0.95 

99100 

2.5 

0.6 

2 

2.30 

1.19 

36530 

0.476 

0.00 

-1.02 

92670 

2.5 

1  0.25 

2.95 

1.26 

30317 

0.403 

0.00 

-0.88 

107070 

2.5 

1 

1 

2.27 

1.33 

35160 

0.531 

0.00 

-0.93 

101800 

2.5 

1 

2 

2.18 

1.36 

35720 

0.558 

0.00 

-0.99 

95830 

5 

0.2  0.25 

11.85 

1.21 

9058 

0.102 

0.00 

-1.43 

65975 

5 

0.2 

1 

3.31 

0.69 

31060 

0.205 

0.65 

-0.65 

93880 

0.785 

5 

0.2 

2 

2.49 

1.33 

33310 

0.492 

0.30 

-0.68 

113100 

1.154 

5 

0.6  0.25 

4.65 

1.17 

21578 

0.247 

0.02 

-1.10 

84777 

1.553 

5 

0.6 

1 

2.60 

1.12 

34350 

0.405 

0.00 

-0.90 

105100 

5 

0.6 

2 

2.36 

1.14 

36400 

0.451 

0.00 

-0.95 

98930 

5 

1  0.25 

2.95 

1.26 

30355 

0.402 

0.00 

-0.85 

111580 

5 

1 

1 

2.32 

1.30 

35170 

0.510 

0.00 

-0.85 

111400 

5 

1 

2 

2.24 

1.32 

35670 

0.531 

0.00 

-0.90 

105400 

*  Mode  Mixity  =  tan'‘(Kn/Ki) 


Table  2.5:  Stress  intensity  factors  in  27  different  joint  configurations  at  the  left  and 
right  crack  tips  for  a  positive  moment  load. 
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Joint 

Left  or 

opening  tip 

*ight  or 

closing  ti 

P 

H/h 

w/h  |  a/h 

* 

fii 

Mc 

Mode 

Mixity* 

 * 

fn 

Mc 

Mode 

Mixity* 

1 

0.2  0.25 

11.84 

1.28 

9053.2 

0.107 

0 

-1.56 

60508 

1 

0.2  1 

2.93 

0.83 

33640 

0.276 

0.745 

-0.831 

76440 

0.840 

1 

0.2  2 

1.79 

0.83 

48830 

0.435 

0.216 

-0.828 

101000 

1.316 

1 

0.6  0.25 

4.54 

1.32 

33810 

0.282 

0 

-1.583 

59620 

•1 

0.6  1 

2.13 

1.30 

36910 

0.547 

0 

-1.298 

72740 

1 

0.6  2 

1.65 

1.40 

40600 

0.703 

0 

-1.397 

67560 

1 

1  0.25 

2.71 

1.40 

31040 

0.477 

0 

-1.587 

59490 

1 

1  1 

1.85 

1.45 

37760 

0.666 

0.071 

-1.45 

63610 

1.522 

1 

1  2 

1.63 

1.45 

40020 

0.727 

0.034 

-1.45 

64300 

1.547 

2.5 

2.5 

2.5 

0.2 

0.2 

0.2 

0.25 

1 

2 

11.89 

2.52 

0.75 

1.14 

0.29 

0.18 

9049 

42430 

134500 

0.095 

0.115 

0.235 

0 

-1.485 

-0.339 

-0.198 

63579 

93170 

102800 

0.319 

0.195 

1.025 

1.005 

2.5 

0.6 

0.25 

4.35 

0.80 

23934 

0.182 

0 

-1.413 

66784 

2.5 

0.6 

1 

1.28 

0.39 

75950 

0.294 

0.348 

-0.586 

123000 

1.035 

2.5 

0.6 

2 

0.17 

0.26 

274000 

0.983 

0.841 

-0.33 

109000 

0.374 

2.5 

1 

0.25 

2.22 

0.62 

44450 

0.274 

0 

-1.41 

66947 

2.5 

1 

1 

0.75 

0.39 

111800 

0.483 

0.328 

-0.736 

104600 

1.152 

2.5 

1 

2 

0.05 

0.27 

321200 

1.376 

0.841 

-0.385 

104200 

0.429 

5 

5 

5 

0.2 

0.2 

0.2 

0.25 

1 

2 

11.88 

2.53 

0.72 

1.14 

0.26 

0.11 

9050.5 

42460 

147100 

0.095 

0.104 

0.157 

0.089  -1.475 
1.026  -0.315 
1.036  -0.135 

62223 

94540 

102700 

1.511 

0.298 

0.130 

5 

0.6 

0.25 

4.36 

0.78 

23918 

0.177 

0  -1.403 

67279 

5 

0.6 

1 

1.27 

0.32 

79140 

0.246 

0.349  -0.542 

130400 

0.999 

5 

0.6 

2 

0.09 

0.13 

530500 

0.951 

0.887  -0.219 

113400 

0.242 

5 

1 

0.25 

2.22 

0.58 

44861 

0.256 

0  -1.401 

67378 

5 

1 

1 

0.71 

0.29 

128500 

0.386 

0.329  -0.678 

111600 

1.119 

5 

1 

2 

0.07 

0.10 

700700 

0.987 

0.883  -0.239 

112400 

0.264 

*  Mode  Mixity  =  tan'‘(Kn/Ki) 


Table  2.6:  Stress  intensity  factors  in  27  different  joint  configurations  at  the  left  and 
right  crack  tips  for  a  negative  moment  load. 
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Printed  on:  Mon  Jul  30  17:08:50  Eastern  Daylight  Time  2007  ft?  ABAQUfl 


Figure  2.14:  Contour  plot  of  joint  configuration  (H/h  =  2.5,  w/h  =  0.6,  a/h  =  2)  with  M' 

loading  illustrating  hinging  effect. 
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become  concentrated  at  the  hinge  point.  This  leads  to  a  visually  unaffected  left  notch  and 
a  rather  high  tensile  stress  concentration  at  the  right  notch.  The  green  boxes  show  that 
some  of  these  configurations  have  closing  notch  critical  moment  lower  than  the  one  at 
opening  notch.  This  shows  that  these  configurations  are  more  susceptible  to  crack 
propagation  at  the  closing  overlap  side  rather  than  the  opening  side  of  the  joint.  This  is 
not  much  of  a  concern,  though,  because  it  will  be  shown  in  Chapter  3  that  while  cracks 
do  propagate  on  the  closing  side  they  require  much  more  steps  to  reach  failure,  than  what 
is  usually  required  at  the  opening  side. 

2.7  Metamodel  Construction  and  Optimization 

In  order  to  optimize  a  laser  welded  lap  joint,  the  fracture  strength  has  to  be 
evaluated  using  the  finite  element  method  for  a  large  number  of  candidate  designs.  The 
computational  cost  of  performing  hundreds  of  detailed  finite  element  analyses  can  be 
prohibitive.  One  way  to  overcome  this  problem  is  to  generate  an  approximation  to  the 
finite  element  values  at  a  much  lower  cost.  Such  approximations  are  often  called 
metamodels  in  that  they  provide  a  “model  of  the  model”. 

Based  on  the  results  of  a  limited  number  of  finite  element  analyses,  simple  curve 
fits  were  initially  tried  using  the  mathematical  software  MATLAB  to  obtain  the  stress 
intensity  factors  for  arbitrary  joint  configurations.  An  exponential  curve  fit  was  tested 
over  the  data  obtained  from  parametric  study  and  the  least  square  errors  were  minimized 
using  a  single-objective  genetic  algorithm.  The  method  did  not  provide  satisfactory 
results  due  to  the  high  nonlinearity  of  the  data.  This  prompted  the  investigation,  and  the 
ultimate  use,  of  metamodels. 
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Taking  a  simple  mathematical  approach  to  metamodels,  let  the  vector  jc  denote  the 
decision  variables  and  y  =  fix)  denote  the  true  output  obtained  from  a  finite  element 
analysis  of  the  joint  at  those  decision  variables.  The  corresponding  computationally 
efficient  metamodel  approximation  is  denoted  by  y  =  g(jc).  The  closeness  of  the  two 
results  fix)  and  g(x)  shows  the  accuracy  of  a  metamodel.  There  are  several  techniques  for 
approximating  the  function  fix)  by  g(jt)  such  as  response  surface  methodology  (RSM), 
kriging,  support  vector  regression,  neural  networks  and  radial  basis  functions  (RBF)  (Jin 
et  al.,  2001;  Fang  et  al.,  2005). 

During  the  method  selection  process  two  popular  methods,  RSM  and  RBF,  were 
compared.  RSM  is  very  simple  and  efficient,  but  is  not  suitable  for  creating  global 
models  over  the  entire  design  space  for  highly  non-linear  problems.  In  a  work  by  Yang  et 
al.  (2002)  it  is  suggested  that  RSM  can  be  ineffective  in  multi-objective  optimization 
problems  because  the  response  region  of  interest  would  never  reduce  to  a  small  area  that 
is  good  for  all  conflicting  objectives.  Jin  et  al.  (2001)  compared  several  metamodels 
using  fourteen  different  problems.  Among  those,  the  RBF  method  was  found  to  be  the 
best  for  both  large  and  small  scale  problems.  In  addition,  the  method  was  also  found  to  be 
best  for  overall  performance  with  respect  to  accuracy,  robustness,  problem  type 
flexibility,  sample  size,  efficiency  and  simplicity.  Krishnamurthy  (2003)  compared  other 
metamodels  with  an  augmented  RBF  and  found  it  to  be  very  promising  and  efficient. 
RBFs,  on  the  other  hand,  are  computationally  more  expensive  because  it  uses  a  series  of 
complex  functions  for  a  single  model.  This  becomes  apparent  when  solving  a  multi¬ 
objective  optimization  problem  where  a  set  of  trade-off  solutions  needs  to  be  obtained. 
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The  RBF  procedure  uses  a  linear  combination  of  basis  functions  centered  at  each 
sampling  point.  This  can  be  represented  as 

g(x)  =  X^(|x-x,||),  (20) 

»=i 

where  n  is  the  number  of  sampling  points,  x  is  the  vector  of  design  variables  ( H/h ,  w/h 
and  a/h  for  the  problem  here),  jc,  is  the  vector  of  design  variables  at  the  /,h  sampling  point, 
<f>  is  the  basis  function  and  X,  are  the  unknown  weighting  coefficients.  Some  of  the  most 
commonly  used  basis  functions  are  the  thin  plate  spline,  Gaussian,  multiquadric,  inverse 
multiquadric  and  logistic  basis  functions.  The  data  generated  in  the  parametric  study  was 
found  to  be  best  suited  for  use  with  the  multiquadric  basis  function  which  is  defined  as 

<f>{r)  =  yfr2  +1.  (21) 

An  augmented  radial  basis  function  method  is  utilized  to  construct  the  required 
metamodel.  The  RBF  defined  in  equation  (20)  did  not  provide  a  satisfactory  model  when 
applied  to  the  sampling  data  and  thus  it  was  augmented  using  exponential  functions  to 
better  approximate  the  global  response  of  the  true  function.  The  augmented  RBF  model  is 
expressed  in  the  form 

S(x)  =  2>(|lx  "  XJ)  +  2/,/?,(x),  (22) 

«=>  y=i 

where  pj  are  the  global  approximation  functions,  c,  are  the  corresponding  unknown 
weight  coefficients  and  m  is  the  number  of  linearly  independent  global  approximation 
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functions  used  to  represent  the  global  response.  The  chosen  global  approximation 
functions  is 

P  =  p  g-¥l  g-«2*2  e~a3X3  e~alX\~a2X2  e~<*2x2-<*3x3  e~<*3X3~«\X\  g -<*2*2  "«3*3  J  T 

where  x,  represent  the  variables  H/h,  w/h  and  a/h  for  /  =  1,2  and  3  respectively  and  a,  are 
corresponding  coefficients  which  are  optimized  so  as  to  provide  as  best  fit  to  the  data  as 
possible.  In  general,  the  stress  intensity  factors  obtained  from  the  finite  element  analyses 
were  found  to  be  smoothly  increasing  or  decreasing  monotonic  functions  of  the  design 
variables.  Therefore,  the  global  response  of  RBF  is  appended  using  exponentially 
increasing  or  decreasing  functions.  By  combining  an  appropriate  global  model  and  the 
RBF  methodology,  a  quality  approximation  is  created  which  possesses  the  desirable  trait 
of  returning  the  true  function  values  at  the  sampling  points.  As  for  the  undetermined 
coefficients,  the  exponents  a*  and  the  weight  coefficients  Cj  of  the  global  approximation 
functions  are  obtained  by  the  method  of  least  squares.  The  sum  of  the  square  of  the  errors 
at  the  sampling  points  is  minimized  using  a  real-coded  single-objective  genetic  algorithm. 
These  three  sets  of  a  and  c  coefficients  will  be  different  depending  on  the  type  of  load, 
the  location  of  the  notch  and  the  fracture  mode.  These  coefficients  will  be  different  if 
only  a  two-dimensional  curve  fit  is  desired  (with  one  design  variable  held  fixed)  instead 
of  the  full  three-dimensional  model.  While  it  may  be  a  concern  for  other  methods,  there  is 
no  need  to  impose  any  orthogonality  condition  to  obtain  c  coefficients.  Once  the  a  and  c 
coefficients  are  obtained  from  the  GA,  the  only  unknown  variables  to  be  determined  are 
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the  weighting  coefficients  A  which  are  determined  by  solving  the  linear  system  of 
equations  (Krishnamurthi,  2003) 

A  A  +  Pc  =  g,  (23) 

where  AtJ=  Ml*r*,||),  Pij  =  Pj(x,)  (i  =  1,2 =  1,2,. ..m),  A  =  [Ah  A2,...,  4]T,  c  =  [cu 
C2,...,  cm]T  and  g  =  [g(*i),  g(x2),...,  g(Xn)]T-  As  previously  alluded  to,  it  should  be 
reiterated  that  the  augmented  RBF  approximation  passes  through  all  the  sampling  points 
exactly.  In  other  words,  the  output  of  the  metamodel  is  identical  to  the  true  function  value 
at  the  sampling  points  unlike  a  least  squares  approximation. 

The  outlined  method  is  used  to  construct  a  metamodel  for  the  stress  intensity 
factor  values  obtained  for  positive  and  negative  moment  loading  cases.  Using  the  data  of 
twenty  seven  sampling  points  for  each  separate  load  case,  the  genetic  algorithm  is  utilized 
to  obtain  the  a  and  c  coefficients.  For  full  three-dimensional  (all  three  design  variables) 
validations,  four  sets  of  three  a  coefficients  and  eight  c  values  are  evaluated  for  each  load 
type.  One  such  three-dimensional  set  of  coefficients  for  the  M  +  load  case  with  stress 
intensity  factors  at  the  right  notch  for  both  mode-I  and  II  is  given  in  Table  2.7.  If  one  of 
the  variables  is  held  constant  and  a  two-dimensional  fit  is  desired  for  further  validation, 
twelve  sets  of  three  a  coefficients  and  four  c  values  are  obtained  for  each  load  type.  One 
of  these  sets  for  a  use  in  a  two-dimensional  evaluation  for  representing  the  stress 
intensities  for  both  fracture  modes  at  the  left  notch  under  M'  load  is  given  in  Table  2.8. 
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No. 

c 

a 

Mode-1 

Mode-ll 

Mode-I 

Mode-ll 

1 

2.2166 

1.32611 

4.15975 

5.21889 

2 

-34.2 

35.3324 

3.40316 

3.10918 

3 

0.2898 

-0.8773 

2.96254 

9.48687 

4 

0.1049 

-1.1726 

5 

-6.358 

-68.995 

6 

42.615 

-10.06 

7 

39.03 

9.28647 

8 

77.755 

200 

Table  2.7:  GA  optimized  c  and  or  coefficients  for  right  notch  stress  intensity  factors 
for  a  full  3D  validation  case  under  M  +  load. 


H/h 

mode 

i 

2 

3 

4 

1 

1 

c 

1.686 

0.028 

0.774 

38.882 

a 

- 

3.423 

2.794 

II 

c 

1.464 

-1.555 

-0.486 

7.049 

a 

_ 

4.452 

6.715 

2.5 

1 

c 

0.156 

0.498 

2.143 

33.018 

a 

- 

3.1294 

2.192 

II 

c 

0.3519 

-0.297 

0.235 

2.220 

a 

- 

1.8538 

2.311 

5 

1 

c 

0.1268 

0.4276 

2.055 

32.571 

a 

- 

3.0642 

2.154 

II 

c 

0.183 

-0.182 

0.261 

1.981 

a 

- 

1.4501 

1.897 

Table  2.8:  GA  optimized  c  and  a  coefficients  for  left  notch  stress  intensity  factors  for 

a  full  2D  validation  case  under  M'  load. 
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2.8  Validation 


As  stated  earlier,  the  augmented  RBF  passes  through  all  sampling  points  exactly 
and  thus  the  only  concern  with  the  RBF  is  whether  or  not  the  model  is  accurate  at 
locations  away  from  the  sampling  points.  Therefore,  the  accuracy  of  the  metamodel  needs 
to  be  checked  at  points  which  do  not  coincide  with  the  sampling  points.  Different 
validation  designs  are  constructed  for  both  two  and  three-dimensional  curve  fits.  In  a 
two-dimensional  plot,  function  values  are  evaluated  by  fixing  one  of  the  three 
dimensionless  ratios  H/h,  w/h  and  a/h  and  changing  the  other  two,  while  in  a  three- 
dimension  validation,  all  three  ratios  are  varied.  Ten  finite  element  analyses  are 
performed  for  each  load  type  and  for  both  the  two  and  three-dimensional  cases  to  obtain 
the  stress  intensity  factors  for  use  in  validating  the  metamodel.  These  results  are  labeled 
as  validation  points.  To  aid  in  the  validation  process,  two  and  three-dimensional  plots  are 
presented  separately  for  two  different  loading  conditions  which  allow  one  to  visualize  the 
error  in  the  metamodel. 

2.8.1  Two-dimensional  Validation 

For  comparisons  concerning  the  effects  of  load  type,  the  metamodel  input  file  can 
be  changed  from  M '  left  (or  opening)  tip  results  to  M+  right  (or  opening)  tip  results. 
However,  the  M  right  tip  plots  will  come  out  to  be  very  similar  to  the  plots  for  the  M' 
left  tip  because,  if  we  look  at  Tables  2.5  and  2.6  for  an  H/h  ratio  of  1,  it  is  noticed  that 
there  are  negligible  changes  in  the  stress  intensity  values.  Therefore,  instead  of  separating 
the  results  based  solely  on  load  type,  the  two-dimensional  results  have  been  further 
divided  into  left  and  right  notch  plots.  For  two-dimensional  cases  shown  here,  only  the 
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left  notch  will  be  studied  for  two  types  of  loads  while  in  the  next  section  for  the  three- 
dimensional  cases,  the  right  notch  will  be  discussed  for  both  load  types. 

Two  dimensional  surface  plots  of  the  normalized  stress  intensity  factors  fK,~  '  and 
fuM'J  obtained  from  the  metamodel  are  shown  in  Figs.  2.15  (a)  and  (b).  The  thickness 
ratio  Hlh  is  chosen  to  be  equal  to  1  and  the  normalized  stress  intensities  for  M '  at  the  left 
notch  are  plotted  as  a  function  of  alh  and  w/h.  The  plot  shows  an  infinite  number  of 
designs  among  which  ten  are  chosen  to  be  validated  using  actual  finite  element  results. 
The  visible  green  asterisks  are  sampling  points  while  the  yellow  asterisks  are  validation 
points.  The  dark  red  portion  of  the  surface  denotes  maximal  values  of  the  function.  It 
would  be  interesting  to  observe  the  differences  in  the  surface  plot  if  the  load  is  reversed 
(i.e.  M  +)  whilst  keeping  other  factors  constant  ( H/h  =  1,  left  notch).  The  plot 
corresponding  to  this  case  is  shown  in  Figs.  2.15(c)  and  (d)  for  mode-I  and  II 
respectively.  Looking  at  the  mode-I  plot  for  M ",  it  can  be  seen  that  the  maximum  stress 
intensities  are  obtained  when  the  joint  overlap  is  at  a  minimum.  For  the  M  +  case,  it  is 
noted  that  maximum  stress  intensities  occur  for  designs  with  the  small  and  mid-length 
non-welded  and  welded  regions  respectively.  With  these  observations  established, 
consideration  is  given  to  the  other  design  variable.  The  large  errors  in  the  mode-I,  M  +, 
left  tip  plot  are  due  to  the  fact  that  the  actual  value  at  these  locations  is  zero  whilst  the 
model  predicts  a  small,  but  finite  value.  More  sample  analyses  are  required  to  better 
capture  the  closing  phenomenon  which  leads  to  this  behavior. 
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Figure  2.15:  Two-dimensional  surface  plots  of  geometric  configuration  factor  values  at  the  left 
notch  for  H/h= 1:  (a)  mode-I  and  (b)  mode-II  under  M~,  (c)  mode-I  and  (d)  mode-II  under  M  +. 
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The  surface  curvature  clearly  illustrates  the  nonlinearity  in  the  data  and  reinforces 
the  importance  of  using  the  metamodeling  approach.  The  percentage  errors  between  the 
detailed  finite  element  analysis  results  at  the  left  notch  and  the  corresponding  metamodel 
predictions  for  H/h=\  and  both  loading  types  are  given  in  Tables  2.9(a)  and  (b).  It  is 
observed  that  the  metamodel  is  fairly  accurate  and  that  the  error  is  on  the  order  of  one  or 
two  percent.  At  some  locations  higher  errors  can  be  observed,  the  reason  for  which  will 
be  addressed  later.  Mode-I  error  is  not  shown  in  Table  2.9(b)  because  even  though  the 
value  provided  is  not  poor  in  a  practical  sense,  the  error  is  technically  infinite. 

2.8.2  Three-dimensional  Validation 

Here  a  three-dimensional  model  is  used  for  explaining  the  workings  of  the 
metamodel.  In  a  three-dimensional  metamodel,  one  can  obtain  values  of  stress  intensity 
factors  for  a  large  combination  of  the  three  dimensionless  ratios.  The  M +  right  notch  case 
is  considered  here  to  explain  the  metamodeling  process.  In  this  case,  the  function  value 
fiM+r  (or  fui,r)  is  evaluated  to  create  mxnxp  volume  array  of  data  for  plotting.  Slice 
volumetric  plots  for  M  +  mode-I  and  II  are  shown  in  Figs.  2.16(a)  and  (b)  where  darker 
red  portions  represents  a  higher  geometric  configuration  factor  value,  and  conversely, 
darker  blue  regions  indicate  a  lower  factor  value.  For  the  purposes  of  this  discussion, 
recall  that  for  the  joint  geometry,  a/h  and  w/h  combined  represent  the  overlapped  area, 
while  H/h  represents  the  lower  plate  thickness.  As  apparent  from  Fig  2.16(a),  the  smaller 
the  total  overlap  area  the  greater  the  tensile  stress  concentration.  This  effectively  reduces 
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(a) 


Test  Configurations 

Mode- 

Mode-ll 

H/h 

w/h 

a/h 

fl” 

%Error 

fit” 

*fn” 

%Error 

1 

0.5 

1.2 

1.999 

2.030 

-1.551 

1.222 

1.262 

-3.241 

1 

0.7 

0.5 

2.844 

2.781 

2.233 

1.349 

1.340 

0.667 

1 

0.8 

1.5 

1.744 

1.750 

-0.367 

1.398 

1.423 

-1.788 

1 

0.9 

1.6 

1.698 

1.700 

-0.118 

1.428 

1.447 

-1.310 

(b) 


Test  Configurations 

Mode-1 

Mode-ll 

H/h 

w/h 

a/h 

ft*4 

f,rj 

*  hr1 

%Error 

1 

0.5 

1.2 

0 

0.0638 

1.224 

1.295 

5.83 

1 

0.7 

0.5 

0 

0.0284 

1.407 

1.407 

0.00 

1 

0.8 

1.5 

0 

0.0571 

1.398 

1.381 

1.22 

1 

0.9 

1.6 

0 

0.0995 

1.43 

1.421 

0.44 

*  Metamodel  results 

Table  2.9:  Two-dimensional  validation  results  checked  at  the  left  notch  for  H/h= 1  for:  M~  and 

(b)  M+  load. 


78 


Figure  2.16:  Three-dimensional  volumetric  slice  plot  for  M+,  right  notch  tip:  (a)mode-I, 

(b)mode-II. 
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the  strength  of  the  joint,  which  conceptually  makes  sense.  For  mode-II,  the  plot  shows 
that  stress  intensities  are  higher  for  geometries  with  high  non-welded  lengths,  which 
represents  the  contact  area  between  the  left  surfaces.  As  this  contact  area  increases,  the 
shear  stress  increases.  If  the  figure  is  rotated  and  one  inspects  the  slices  on  the  vertical 
axis,  it  is  observed  that  the  shear  stress  distribution  is  highest  when  the  lower  plate 
thickness  is  at  a  minimum.  This  is  because  large  lower  plate  thicknesses  possess  a  high 
stiffness  which  reduces  the  shear  stress  concentration. 

Table  2.10  gives  the  percentage  errors  between  the  finite  element  analyses  for  ten 
different  joint  configurations  and  the  corresponding  metamodel  results  for  the  M  '  right 
notch  case.  It  is  seen  that  the  errors  are  within  seven  percent  for  all  validation  points, 
which  for  such  a  small  number  of  sampling  points,  is  acceptable. 

As  seen  in  Table  2.6,  there  are  some  joint  configurations  which  show  low  failure 
strength  on  the  closing  side  of  the  overlap.  This  ascertains  the  requirement  of  forming  a 
metamodel  for  the  closing  side  also.  Data  has  been  extracted  from  the  same  twenty  seven 
analyses  and  ten  validation  points  to  construct  data  for  the  closing  notch.  There  are  two 
load  types  under  consideration.  For  M  +,  the  left  edge  of  the  weld  is  the  closing  notch 
while  for  M ',  the  right  notch  will  close.  Tables  2.1 1  (a)  and  (b)  present  the  data  for  a 
two-dimensional  metamodel  validation  for  the  left  notch  with  M  4  and  for  a  three- 
dimensional  validation  of  the  right  notch  with  M  \  respectively.  The  errors  have  not  been 
shown  due  to  the  often  indeterminate  nature  of  the  error  for  mode-I.  However,  it  can  be 
seen  that  due  to  the  very  low  values  of  the  geometric  configuration  factors,  any  slight 
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Test  Configurations 

Mode-1 

Mode-ll 

H/h 

w/h 

a/h 

/r,r 

| 

%Error 

.fuM+'r 

\f,rr 

%Error 

1.2 

0.5 

1.2 

2.20 

2.3261 

-5.612 

1.20 

1.1625 

2.801 

1.5 

0.7 

0.5 

3.09 

3.0281 

2.003 

1.23 

1.1855 

3.303 

1.8 

0.8 

1.5 

2.19 

2.27 

-3.700 

1.33 

1.2614 

5.015 

2 

0.9 

1.6 

2.18 

2.2326 

-2.554 

1.35 

1.2951 

4.350 

2 

0.4 

0.4 

5.14 

5.3332 

-3.739 

0.99 

1.0628 

-6.922 

2.2 

0.5 

1.1 

2.57 

2.5863 

-0.713 

1.10 

1.0802 

1.442 

2.4 

0.9 

0.6 

2.61 

2.5462 

2.283 

1.26 

1.2815 

-1.787 

3 

0.7 

0.5 

3.16 

3.0602 

3.112 

1.15 

1.1726 

-2.321 

3.5 

0.4 

0.5 

4.42 

4.5666 

-3.429 

0.94 

0.9675 

-2.926 

4.2 

0.8 

0.6 

2.77 

2.6664 

3.723 

1.19 

1.1219 

6.039 

*  Metamodel  results 


Table  2.10:  Comparison  of  Metamodel  with  finite  element  results  with  Ten  three- 
dimensional  validation  designs  for  the  M  +  right  notch  case  for  both  mode-I  and  II. 
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change,  even  in  the  second  decimal  place,  makes  a  big  difference  in  percentage  error.  In 
some  analysis  types  for  both  the  two  and  three-dimensional  spaces,  comparisons  of  the 
validation  points  with  the  metamodel  values  give  high  errors.  These  results  have  been 
checked  for  different  basis  functions  of  the  metamodel  as  well,  and  it  has  been  found  that 
multiquadric  works  the  best.  Therefore,  the  above  stated  errors  could  potentially  be  due 
to  changes  in  the  mesh  partition  size  used  for  each  of  the  validation  designs.  As  has  been 
previously  discussed,  changes  in  the  partition  size  have  a  significant  impact  on  the 
obtained  results.  Another  potential,  and  probable,  cause  could  be  the  small  number  of 
sampling  points  used  in  constructing  the  metamodel.  With  the  nonlinearity  present  in  the 
data,  more  sampling  points  would  definitely  improve  the  results  albeit  at  the  cost  of  a 
greatly  increased  workload  in  terms  of  analysis  runs  and  post-processing. 

Although  they  are  a  little  tedious  to  construct  and  give  moderate  errors  for  some 
configurations,  metamodels  provide  a  means  of  closely  approximating  the  stress  intensity 
factors  for  arbitrary  joint  configurations  in  a  fraction  of  the  time  that  it  takes  to  perform  a 
detailed  finite  element  analysis. 
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(a) 


Test  Configurations 

Mode-I 

Mode-ll 

H/h  | 

w/h 

a/h 

f,M+’r  1 

*frr 

fnm'  1 

*fuM+r 

1 

0.5 

1.2 

0 

0.0638 

1.224 

1.295 

1 

0.7 

0.5 

0 

0.0284 

1.407 

1.407 

1 

0.8 

1.5 

0 

0.0571 

1.398 

1.381 

1 

0.9 

1.6 

0 

0.0995 

1.43 

1.421 

2.5 

0.4 

0.4 

0 

0.0738 

0.935 

1.049 

2.5 

0.5 

1.1 

0 

0.0389 

0.937 

0.924 

2.5 

0.9 

0.6 

0 

-0.0052 

0.876 

0.964 

2.5 

0.7 

0.5 

0 

0.0073 

0.894 

1.048 

5 

0.4 

0.5 

0 

0.1256 

0.837 

0.898 

5 

0.8 

0.6 

0 

-0.01 

0.837 

0.928 

(b) 


Test  Configurations 

Mode-I 

Mode-ll 

H/h 

w/h  [ 

a/h 

/r,r  | 

•k  fM+J 

fnM+r 

* 

1.2 

0.5 

1.2 

0.23 

0.3034 

0.94 

1.012 

1.5 

0.7 

0.5 

0 

0.1807 

1.114 

1.212 

1.8 

0.8 

1.5 

0.56 

0.6064 

0.643 

0.733 

2 

0.9 

1.6 

0.64 

0.6701 

0.58 

0.634 

2 

0.4 

0.4 

0 

0.1377 

1.027 

1.167 

2.2 

0.5 

1.1 

0.454 

0.4961 

0.525 

0.537 

2.4 

0.9 

0.6 

0 

0.2042 

0.996 

0.951 

3 

0.7 

0.5 

0 

0.1846 

1.003 

0.983 

3.5 

0.4 

0.5 

0 

0.2345 

0.846 

0.957 

4.2 

0.8 

0.6 

0 

0.2947 

0.932 

0.896 

*  Metamodel  results 


Table  2.11:  Metamodel  validation  results  checked  at  ten  joint  configurations  for:  (a)  2D 
cases  at  left/closing  notch  for  both  mode-I  and  II,  and  (b)  3D  M '  cases  at  right/closing 

notch  for  both  mode-I  and  II. 
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2.9  Results  and  Discussion 


In  this  section,  optimal  design  studies  are  performed  on  two  model  problems  and 
the  results  are  presented.  These  design  studies  will  demonstrate  the  procedure  of  using 
the  established  model.  The  joint  geometry  is  similar  to  the  lap  joint  shown  in  Fig  2.1.  For 
both  studies,  the  insert  and  sub-frame  are  assumed  to  be  made  of  HSLA-65  steel.  The 
yield  strength  and  critical  stress  intensity  factors  are  taken  to  be  Sy  =  65  Ksi  and  K|C  = 
1 13.76  ksiVin.  Both  of  these  model  problems  utilize  previously  constructed  metamodels. 


2.9.1  Design  Problem  I:  Optimal  Weld-Width  with  a  Positive  Moment 

We  consider  the  design  of  a  laser  welded  joint  that  is  subjected  to  a  positive 
moment  of  Mdesign  =  750  lb-in/in.  The  stress  in  the  insert  material  is  obtained  from  beam 
bending  theory: 


M 


6M 
h  2  ' 


(24) 


The  factor  of  safety  against  yielding  5y/a*  is  chosen  to  be  2.  Obtaining  a*  from  this  yield 
factor  of  safety  expression  and  plugging  the  result  into  equation  (24)  gives  the  thickness 
of  the  insert  material.  Rounding  the  obtained  thickness  to  the  next  standard  dimension 
results  in  a  value  of  h  =  3/8".  Assuming  a  sub-frame  thickness  of  H  =  1 ",  the  ratio  of  the 
insert  thickness  to  the  sub-frame  thickness  is  H/h  =  2.67.  The  metamodel  for  the  right 
notch  positive  moment  case  is  employed  to  compute  the  geometric  configuration  factors 
flM+,r(2.67,f,f)  for  mode-I  and  II.  If  the  ratio  a/h  is  fixed  at  2,  a  range  of  values 

varying  with  only  w/h  could  be  obtained.  The  stress  intensity  factors  K^+'r  and  KuM+  r  at 
the  right  notch  are  evaluated  using  equation  (18)  where  M  is  replaced  by  Mdesign.  The 
crack  initiation  angle  0o  is  computed  by  inserting  the  obtained  K/  and  Kn  values  in 
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equation  (6).  With  90,  Ki  and  Ku,  the  critical  fracture  moment  act  is  obtained  using  the 
maximum  tensile  stress  criteria  for  mixed  mode  failure  via  equation  (8).  Once  the  critical 
moment  for  the  model  is  obtained,  the  factor  of  safety  against  fracture  can  be  defined  as 
iYfr  =  Mfract/Mdesign-  Fig.  2.17  shows  the  fracture  safety  factor  as  a  function  of  the  weld- 
width  ratio  w/h  for  several  different  choices  of  a/h.  It  is  noticed  that  the  factor  of  safety  is 
greater  than  2  for  the  range  of  joint  dimensions  considered  in  this  study.  When  the  width 
of  the  non-welded  region  a/h  <=  1.0,  the  factor  of  safety  increases  as  the  width  of  the 
non-welded  region  w  increases.  The  trend  reverses  for  values  of  a/h  >  1.0.  For  a  given 
weld-width  w,  it  is  found  that  a  significant  increase  in  the  fracture  safety  factor  can  be 
achieved  for  an  a/h  value  of  approximately  1.5.  If  a/h  is  larger  than  1.5,  the  weight  of  the 
joint  would  increase  whereas  the  strength  would  remain  virtually  unchanged  or  decrease 
slightly.  This  drop  in  the  strength  value  can  be  explained  by  the  mode  mixity  factor, 
shown  in  Figure  2.18.  If  mode  mixity  is  plotted  versus  w/h,  it  can  be  seen  that  its  value 
increases  for  every  increment  in  the  a/h  ratio.  From  this  one  can  infer  that  the  mode-11 
factor  plays  a  more  significant  role  in  failure  for  higher  overlaps.  As  the  overlap 
continues  to  increase  after  a  certain  value  of  a/h,  more  shear  stresses  develop  which  lead 
to  lowering  the  joint  strength  and  fastening  the  failure  of  the  joint. 
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Safe  designs  for  different  weld-widths 


Figure  2.17:  Fracture  safety  factor  as  a  function  of  weld-width  for  different  overlaps  with  M  +, 

right  tip  function  evaluations. 

Influence  of  K^/  on  failure 


Figure  2.18:  Mode  mixity  as  a  function  of  weld-width  for  different  overlaps  with  M  +,  right  tip 

function  evaluations. 
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Using  these  results,  the  optimum  joint  dimensions  for  this  design  case  study  are  h 
=  0.375  in,  w  =  (0.2  x  h  )  =  0.075  in  and  a  =  (1.5  *  h)  =  0.5625  in.  This  model  problem 
demonstrates  the  benefit  of  performing  a  limited  number  of  detailed  finite  element 
simulations  to  aid  in  the  construction  of  a  metamodel  that  can  then  be  later  used  to 
optimize  laser  welded  lap  joints. 

2.9.2  Design  Problem  II:  Optimal  Weld-Width  with  a  Negative  Moment 

As  was  discussed  earlier,  there  are  some  non-zero  mode-I  stress  intensity  factors 
at  the  closing  side  of  the  joint,  most  notably  in  the  negative  moment  cases.  Therefore,  for 
this  design  problem,  a  negative  moment  is  applied  and  the  critical  load  of  failure  at  the 
left  and  right  notches  of  the  joint  are  analyzed  and  compared.  The  joint  dimensions  are 
identical  to  those  used  in  design  problem  1 .  The  applied  moment  for  this  case  is  Mdesign  - 
-750  lb-in/in.  As  in  the  previous  design  study,  keeping  the  yield  factor  of  safety  at  2  gives 
a  thickness  of  h  =  3/8"  and  an  H/h  value  of  2.67. 

As  described  in  the  previous  design  problem,  the  metamodel  is  utilized  to 
compute  the  factor  of  safety  against  fracture  N&  =  MfraC, /Assign  for  both  the  left  and  right 
notches  when  the  joint  is  subjected  to  a  negative  moment.  These  two  values  are  compared 
and  the  minimum  is  selected  as  the  final  value  for  the  fracture  factor  of  safety.  To  get  a 
better  handle  on  the  problem,  the  range  of  a/h  is  subdivided  into  very  small  segments  and 
a  three-dimensional  surface  plot  is  shown  in  Figs.  2.19(a)  and  (b). 

From  these  plots  it  can  be  observed  that  at  a  minimum  weld-width  of  w/h  =  0.2 
the  factor  of  safety  increases  as  the  overlap  increase  up  to  a  value  of  a/h  =1.5,  beyond 
which  the  value  remains  relatively  unchanged.  Once  the  weld  length  increases  to  0.3h,  a 
slow  increase  in  NfT  is  seen  even  beyond  a/h  =  1.5.  Once  the  curve  peaks,  the  value  of  N(r 
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begins  to  diminish  with  increasing  weld-widths.  A  very  prominent  feature  of  these  plots 
is  the  line  of  peak  value.  This  peak  line  remains  almost  unchanged  for  an  entire  series  of 
configurations,  which  is  rather  interesting.  These  values  possess  total  overlap  ratios  ( a/h 
+  w/h)  in  the  range  of  1.8  to  2.1.  However,  the  maximum  strength  is  seen  away  from  this 
line  at  design  values  of  w/h  =  0.3  and  a/h  =  2.0. 

Although  there  are  a  range  of  designs  to  choose  from,  if  the  errors  discussed  in 
evaluating  the  stress  intensity  factors  in  Section  2.9  are  taken  into  consideration,  and  if 
one  wants  to  check  the  design  with  maximum  possible  factor  of  safety,  one  would  be 
better  off  with  an  a/h  ratio  of  1.65  and  w/h  value  of  0.3.  Thus,  the  optimum  joint 
dimensions  for  this  design  case  study  are  taken  as  h  =  0.375  in,  w  =  (0.3  x  h  )  =  0.1 125  in 
and  a  =  (1.65  *  h  )  =  0.6187  in.  The  factor  of  safety  corresponding  to  this  joint  is  Nfr  = 
39.29.  With  this  design  study  we  saw  that,  even  for  a  complicated  loading  condition  on  a 
lap  joint,  metamodels  prove  to  be  a  very  useful  and  computationally  inexpensive  tool  for 
designing  optimal  lap  joint  configurations. 
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safe  oes  grs  (or  afferent  #ed-*KJtns 


Figure  2.19:  Two-dimensional  surface  plots  of  minimum  fracture  factor  of  safety  at  H/h= 2.67 
compared  between  the  left  and  right  notch  forM-  cases. 
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3.  FATIGUE  CRACK  PROPAGATION 


3.1  Introduction 

In  this  section,  the  two-dimensional  crack  propagation  simulation  software 
FRANC2D  is  used  to  study  crack  growth  in  laser  welded  lap  joints.  Sequential  crack 
propagation  analysis  is  performed  to  obtain  the  crack  trajectories  and  corresponding 
stress  intensity  factors  for  lap  joint  configurations  under  three  different  loads;  positive 
moment,  negative  moment  and  axial  tension.  The  results  of  these  analyses  for  each  load 
case  are  utilized  in  conjunction  with  the  Paris  law  to  obtain  a  fatigue  response  estimate  of 
the  joints.  Two  joint  configurations  are  selected  to  be  studied  under  each  load  case 
among  which  one  joint  is  set  to  be  common  to  all  the  load  cases. 

3.2  Modeling  in  FRANC2D 

Since  the  FEA  program  ABAQUS  v6.6  does  not  provide  a  means  for  modeling 
crack  propagation,  a  fracture-specific  software  is  selected  for  use  in  this  chapter. 
FRANC2D  (FRacture  ANalysis  Code  for  2D  geometries)  is  a  two-dimensional  finite 
element  based  program  for  simulating  curvilinear  crack  propagation  in  plane  stress,  plane 
strain  and  axisymmetric  structures.  The  software  was  developed  by  Professor  Ingraffea 
and  his  graduate  students  at  Cornell  University.  FRANC2D  requires  an  initial  mesh  of 
the  geometry  which  has  to  be  constructed  separately  since  it  does  not  feature  a  mesh 
generation  module.  Although  it  can  be  used  with  input  files  made  in  several  different 
finite  element  software  packages,  the  recommended  pre-processing  software  CASCA  is 
used  to  define  the  joint  geometries  and  generate  the  initial  mesh  for  all  studies  conducted 
in  this  chapter.  Translation  codes  are  also  available  on  the  Cornell  Fracture  Group 
website  to  convert  other  commercial  finite  element  pre-processed  data  into  a  format  that 
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is  suitable  for  FRANC2D  input.  Several  different  mesh  partition  strategies  were 
examined,  as  will  be  discussed  in  Section  3.3.1,  and  the  selected  mesh  strategy  for  one 
such  joint  configuration  is  shown  in  Fig.  3.1. 

3.2.1  Preprocessing 

To  begin  the  process,  a  FE  model  is  created  without  cracks  using  the  external 
mesh  generator  CASCA,  as  shown  in  Fig.  3.1(b).  After  loading  the  CASCA  generated 
file,  the  remaining  pre-processing  is  performed  within  FRANC2D.  Similar  boundary 
conditions  are  defined  in  FRANC2D  as  were  in  ABAQUS  models.  To  simulate  the 
boundary  conditions  in  FRANC2D,  translational  degrees  of  freedom  in  x  and  y  directions 
at  the  left  edge  of  the  lower  plate  are  set  to  be  equal  to  zero.  Point  and  edge  loads  have 
been  applied  at  the  free  edge  of  the  upper  plate  to  create  the  required  couple  and  tensile 
forces,  respectively.  The  next  step  is  to  define  the  material  property  of  the  joint.  The 
entire  model  is  defined  to  be  made  of  any  type  of  steel  with  a  Young’s  Modulus  of  E  = 
29.7x1 03  ksi  and  a  Poisson’s  ratio  of  v  =  0.29.  In  FRANC2D,  contact  modeling  requires 
the  use  of  gap  elements  and  thus  a  second  material  property  set  for  gap  elements  is 
required.  The  gap  element  property  definition  lets  the  user  specify  the  surface  friction 
and  the  interpenetration  allowed  between  the  surfaces  in  contact.  These  elements  are 
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(a) 


Welded  region 


(b) 


Figure  3.1  Mesh  partitioning  scheme  on  joint  configuration  (H/h= 1,  h»//i=0.6,  a/h—  1)  used 
for  almost  all  joint  configurations  analyzed  with  FRACN2D. 
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allocated  at  the  crack  surfaces  to  enforce  contact  conditions.  For  all  the  joints  under 
consideration,  the  problem  type  is  selected  to  be  plane  strain.  Although  there  are  several 
other  options  in  the  program,  only  those  pertaining  to  the  work  conducted  in  this  chapter 
are  mentioned  here.  Once  the  pre-processing  is  done,  one  has  to  introduce  the  non- 
welded  lengths  as  cracks. 

3.2.2  Creation  of  cracks 

In  FRANC2D,  the  user  can  change  the  geometry  as  well  as  the  mesh  of  the 
structure  which  is  essential  for  performing  the  desired  fracture  analyses.  It  may  include 
addition  or  deletion  of  elements,  introduction  or  propagation  of  the  crack,  and  creation  of 
material  interfaces  and  contact  surfaces.  For  the  current  purpose,  two  new  cracks  (to 
model  the  non-welded  overlapped  lengths)  are  introduced  in  the  model  on  either  side  of 
the  known  welded  length.  Cracks  can  be  non-cohesive  (traction  free  crack  surfaces), 
cohesive  or  those  with  gap  elements.  Gap  elements  are  utilized  in  the  present  analysis  to 
model  contact  between  the  plates  in  the  non-welded  overlap  region.  The  entire  process  of 
crack  creation  is  graphically  shown  in  Fig.  3.2  (a),  (b),  (c),  (d)  and  (e).  Upon 
specification  of  an  edge  crack,  the  program  deletes  the  elements  around  the  crack  path, 
creates  crack  surfaces  with  gap  elements,  forms  a  rosette  of  eight  six-node  triangular  T6 
elements  around  the  crack-tip,  and  later  meshes  the  deleted  region  using  the  quad-tree 
algorithm,  which  produces  a  good  mesh  consisting  of  T6  elements  around  the  crack. 
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(e) 

Figure  3.2  Introduction  of  an  edge  crack  in  the  joint  (H/h  =  1,  w/h  =  0.6,  a/h  =  1)  where, 
(a)section  of  Fig.3.1(b),  (b)deletion  of  surrounding  elements,  (c)creation  of  crack  surfaces, 
(d)creation  of  the  rosette  at  the  crack  tip,  (e)complete  mesh. 
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3.2.3  Analysis 

Two  types  of  solvers  are  available  in  FRANC2D  namely,  direct  stiffness  and 
dynamic  relaxation.  Direct  stiffness  is  a  standard  Gauss  elimination  solver  for  linear 
problems,  while  dynamic  relaxation  is  used  for  problems  containing  gap  elements  (or 
rather  any  geometric  nonlinearity,  such  as  contacts).  A  convergence  tolerance  has  to  be 
specified  for  the  dynamic  relaxation  algorithm.  Smaller  tolerances  generally  yield  more 
accurate  results  at  the  expense  of  computational  time. 

3.3  Convergence 

To  ensure  that  the  mesh  is  sufficiently  refined  and  that  the  results  have  converged, 
the  center  crack  model  previously  analyzed  using  the  program  ABAQUS  is  made  and 
meshed  in  CASCA.  Two  types  of  element  distributions  are  seen  while  creating  mesh  in 
CASCA  and  FRANC2D,  namely,  number  of  elements  around  the  crack  (EAC)  and  the 
number  of  elements  over  the  crack  length  (EOCL).  EAC  is  defined  as  the  number  of 
element  seeds,  at  one  edge  of  the  innermost  rectangular  partition  of  the  CASCA  mesh, 
around  the  crack  tip  and  EOCL  defines  the  number  of  elements  along  the  crack  length 
(during  mesh  modification  in  FRANC2D),  which  is  shown  in  Fig.  3.2(d).  As  the  effect 
of  these  element  density  parameters  on  the  final  analysis  result  is  not  known  a  priori,  the 
following  studies  are  conducted;  EOCL  is  fixed  and  EAC  is  increased,  EAC  is  fixed  and 
EOCL  is  increased,  and  both  EAC  and  EOCL  are  increased  simultaneously.  The  results 
are  shown  in  Table  3.1(a).  When  EOCL  is  fixed  at  40  and  EAC  is  increased  from  20  to 
50  elements,  the  difference  between  the  obtained  stress  intensity  factors  and  analytical 
results  remained  almost  unchanged.  Similar  results  are  obtained  in  the  third 
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(a) 


EOCL 

EAC 

MCC'  K, 

Analytical  K, 

Difference  (%) 

40 

20 

77.79 

79.27 

1.867 

40 

30 

77.93 

79.27 

1.690 

40 

40 

77.93 

79.27 

1.690 

40 

50 

77.91 

79.27 

1.716 

EOCL 

EAC 

MCC’  K, 

Analytical  Ki 

Difference  (%) 

20 

20 

78.58 

79.27 

0.870 

25 

25 

78.61 

79.27 

0.833 

30 

30 

78.54 

79.27 

0.921 

40 

40 

78.52 

79.27 

0.946 

50 

50 

78.44 

79.27 

1.047 

EOCL 

EAC 

MCC’  K, 

Analytical  K, 

Difference  (%) 

10 

50 

79.3 

79.27 

0.038 

20 

50 

79.13 

79.27 

0.177 

25 

50 

79.08 

79.27 

0.240 

30 

50 

78.97 

79.27 

0.378 

40 

50 

78.66 

79.27 

0.770 

50 

50 

78.44 

79.27 

1.047 

EOCL 

EAC 

MCC’  K, 

Analytical  Ki 

Difference  (%) 

10 

40 

79.03 

79.27 

0.303 

20 

40 

78.98 

79.27 

0.366 

30 

40 

78.76 

79.27 

0.643 

40 

40 

78.52 

79.27 

0.946 

50 

40 

78.24 

79.27 

1.299 

(b) 


*MCC  is  Modified  Crack  Closure  Method  of  computing  SIFs. 
**Unit  of  stress  intensity  factors  listed  in  the  table  is,  psiVin 


Table  3.1(a)  Convergence  study  using  center  crack  model,  (b)  Comparison  of  FRANC2D  and 
ABAQUS  results,  for  a  representative  joint  configuration  (H/h  =  1,  w/h  =  0.6,  a/h  =  1). 
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case  by  increasing  both  the  element  distribution  parameters  from  20  to  50.  For  the 
second  study,  when  EAC  is  fixed  at  50  elements  it  is  observed  that  the  difference 
increases  with  EOCL  being  increased  from  10  to  50  elements.  Thus,  the  best  result 
occurs  at  an  EOCL  of  10.  The  same  procedure  is  repeated  for  EAC  numbers  of  40,  30 
and  20  with  the  EOCL  element  number  increasing  from  10  to  50.  It  is  clearly  seen  that 
with  an  increase  in  the  EOCL  parameter,  the  difference  increases.  Therefore,  for  almost 
all  meshes  generated  in  this  chapter,  an  EOCL  of  10  is  selected. 

Once  an  appropriate  number  of  elements  to  be  used  around  the  crack  tip  and 
elements  over  the  crack  length  are  determined  using  center  crack  FE  model,  the  next  step 
is  to  make  sure  that  these  modeling  parameters  also  provide  accurate  results  in  the  lap 
joint  models.  To  examine  this,  several  mesh  convergence  partitioning  schemes  were 
analyzed  for  a  joint  configuration  of  H/h  =  1,  w/h  =  0.6,  and  a/h  =  1.  The  two 
partitioning  schemes  used  are  displayed  in  Figs.  3.3(a)  and  (b).  The  mesh  partitioning 
scheme  plays  a  significant  role  in  the  numeric  value  of  the  output.  It  is  expected  that  the 
finer  the  element  mesh  the  more  accurate  the  stress  distribution  results  will  be. 

The  model  shown  in  Fig.  3.3(a)  has  an  EAC  of  80  and  EOCL  of  10.  It  exhibits  a 
difference  of  6%  in  the  Kj  value  and  4%  in  the  Kn  value  when  compared  to  those 
obtained  using  ABAQUS  for  the  same  joint  geometry.  The  model  shown  in  Fig.  3.3(b)  is 
partitioned  in  such  a  way  as  to  obtain  a  more  refined  mesh  around  the  crack  tips  with  an 
EAC  of  20  and  EOCL  of  10.  Table  3.1(b)  compares  results  obtained  using  FRANC2D 
for  the  mesh  distribution  shown  in  Fig.  3.3(b)  with  those  obtained  from  ABAQUS  for  a 
joint  of  dimensions  h  =  H=  a  =  1,  w  =  0.6  subjected  to  a  unit  moment.  The  results 


97 


(a) 


Figure  3.3:  Mesh  partitioning  schemes  of  (a)80  elements  and  (b)20  elements  on  one  side  of 
the  innermost  rectangular  section,  on  the  joint  configuration  (H/h  =  1,  w/h  =  0.6,  a/h  =  1). 
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given  in  the  table  demonstrate  that,  for  this  mesh  scheme,  tighter  convergence  tolerances 
for  the  dynamic  relaxation  algorithm  is  really  not  beneficial  and  an  EOCL  of  1 0,  an  EAC 
of  20  with  a  tolerance  of  1><1 0  6  gives  sufficiently  accurate  results.  As  it  provides 
satisfactory  results  and  reduces  computer  run  time,  a  similar  mesh  partitioning  scheme,  as 
the  one  shown  in  Fig.  3.3(b)  with  four  levels  of  partitions,  is  selected  for  the  joint 
configurations  considered  here. 

3.4  Crack  Propagation  Method 

In  Chapter  2,  the  static  fracture  strength  of  the  joint  was  obtained  by  studying 
fracture  initiation  at  the  notches.  In  addition  to  the  static  fracture  strength,  it  is  helpful  to 
understand  crack  growth  and  corresponding  residual  strength  of  the  joint.  Accordingly,  a 
sequential  crack  propagation  analysis  is  performed  using  FRANC2D  in  this  section.  A 
failure  criterion  which  uses  stress  intensity  factors  to  determine  crack  initiation  angle 
under  mixed  mode  conditions  is  discussed  in  Section  2.4.  As  an  extension  to  that 
methodology,  crack  propagation  can  be  considered  as  a  step-by-step  fracture  process 
using  the  maximum  tensile  stress  failure  criterion.  Crack  propagation  analysis  begins 
with  a  previously  defined  crack.  There  are  generally  three  methods  to  propagate  a  crack 
using  FRANC2D,  which  are  as  follows;  the  standard  method  which  allows  the  user  to 
specify  a  crack  length  increment  and  direction,  the  interface  method  which  is  used  when 
the  crack  tip  is  located  at  a  bimaterial  interface,  and  the  automatic  method.  For  the 
purposes  of  this  work,  the  standard  crack  propagation  algorithm  is  chosen  as  it  is  the  most 
appropriate. 

At  each  step  of  the  crack  propagation  process,  the  stress  intensities  are  first 
obtained  by  the  Modified  Crack  Closure  technique.  These  stress  intensity  factor  values 
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are  used  to  calculate  the  increment  in  crack  length  which  is  then  supplied  to  the  program 
manually  (refer  Section  3.5).  Referring  to  the  method  outlined  in  Section  2.4,  FRANC2D 
then  evaluates  the  crack  propagation  angle  for  each  crack  using  the  maximum  tensile 
stress  criterion.  Subsequently,  the  changes  in  crack  geometries  are  incorporated  into  the 
existing  finite  element  model  and  the  elements  around  the  new  cracks  are  regenerated. 
Thereafter,  the  model  is  analyzed  by  the  dynamic  relaxation  algorithm  and  the  new  stress 
intensity  factors  at  the  new  crack  tips  are  probed  from  the  program  output.  This 
procedure  is  repeated  until  a  crack  surface  reaches  very  close  to  a  boundary  or  another 
crack  surface.  Fig.  3.4  shows  the  modified  joint  configuration  during  the  sequential  crack 
propagation  process  at  the  3rd,  8th  and  16th  step  for  the  joint  JC:  1  -0.6-1  under  axial 
tension.  The  entire  stress  intensity  factor  history  which  is  obtained  as  a  function  of  the 
crack  length  is  then  utilized  in  the  fatigue  analysis  of  the  joint  as  described  in  Section  3.6. 
In  this  whole  process  of  crack  growth,  the  propagation  angle  is  independent  of  the 
fracture  toughness  and  yield  strength,  as  also  shown  in  Section  2.4,  Equation  (6). 
Therefore,  there  would  not  be  any  change  in  the  trajectory  for  different  steels  within  the 
range  of  applicability  of  LEFM.  However,  it  should  be  noted  that  the  residual  strength  of 
the  joint  (which  can  also  be  called  as  critical  moment  of  the  joint)  depends  on  the  fracture 
toughness  K/c  of  the  material. 
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(C) 


Figure  3.4  Crack  propagation  process  at  3,  8  and  16  increment  steps  (a  -  c 
respectively)  on  the  joint  JC:l-0.6-l  under  P+  load. 
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Except  for  one  joint  geometry,  which  is  common  among  all  the  load  cases,  an 
additional  joint  configuration  has  also  been  examined  for  each  load  case.  All  reported 
quantities  discussed  later  in  this  chapter  are  normalized  with  respect  to  either  the  upper 
plate  thickness  h,  applied  moment  load  M  or  tensile  load  P,  critical  stress  intensity  factor 
Kjc,  other  material  constants,  or  a  combination  of  these,  and  are  denoted  with  bar  over  the 
quantity.  These  normalized  quantities  are  listed  hereunder, 

•  Crack  extension, 


6  h’ 

Stress  intensity  factor  for  applied  moment  or  tensile  load, 


3 

h2 

K‘  =  K‘M- 


OR 


1 

h.2 

Ki  =  KtT 


•  Critical  moment  or  residual  strength, 


_  Mc 

Mc  = 

h2.K , 


lc 


Critical  tensile  force  or  residual  strength, 


Pc 


pc=  — 
h2.KIC 

where,  subscript  i  denotes  the  fracture  mode.  To  simplify  and  ease  the  understanding  of 


the  results,  the  following  terminology  is  used  to  define  the  joint  geometries. 


Joint  Configuration 

H  w  a 

—  — .Location 
h  h  h 


In  the  box  above,  the  ratios  following  the  colon  are, 

H/h  =  ratio  of  plate  thicknesses, 

w/h  =  ratio  of  the  weld-width  to  the  upper  plate  thickness,  and 
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a/h  =  ratio  of  the  non- welded  overlapped  length  to  the  upper  plate  thickness. 

Here  ‘Location’  is  used  to  indicate  results  for  the  right  or  the  left  notch  using  the  letters 
‘R’  or  ‘L’,  respectively.  For  example,  JC:2.5-0.2-l,L  would  indicate  the  results 
correspond  to  the  left  notch  of  joint  configuration  H/h  =  2.5,  w/h  =  0.2  and  a/h  =  1. 


3.5  Power-Law  Crack  Growth  Model 


In  a  body  with  multiple  cracks,  such  as  the  lap  joint  configuration  discussed  here, 
it  is  assumed  that  each  crack  will  propagate  by  different  amounts  which  would  depend  on 
their  stress  intensity  factors.  The  following  power-law  equation  is  commonly  employed 
(FRANC3D  manual  v2.6)  to  determine  the  change  in  crack  lengths  at  each  step  of  the 
crack  propagation  process, 


(Ke,r) 

(K^r  ■ 


(25) 


where,  i  is  the  crack  number,  Ae  is  the  prescribed  maximum  normalized  crack  increment 
for  the  step,  is  the  effective  stress  intensity  factor  of  the  ith  crack,  K^x  is  the  largest 
value  of  the  effective  stress  intensity  factor  of  all  cracks  and  b  is  a  power  exponent. 
Considering  an  example,  the  law  states  that  if  there  are  two  cracks  and  the  maximum 
effective  stress  intensity  factor  occurs  at  crack  1,  then  crack  1  will  propagate  by  the 
specified  maximum  value  Ae  while  crack  2  will  extend  by  a  fraction  of  the  specified 
increment  as  computed  by  the  preceding  equation.  If  the  power-law  exponent  b  is  chosen 
to  be  a  large  value  then  crack  2  will  propagate  by  a  smaller  amount.  The  power-law 
exponent  b  is  chosen  to  be  3.0  in  the  present  analysis.  The  selection  of  this  value  is 
discussed  in  detail  in  the  fatigue  section  (Section  3.6)  of  the  chapter.  It  should  be  noted 
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that  this  power-law  only  gives  information  about  the  amount  by  which  a  crack  will 
propagate  and  does  not  produce  the  stress  intensity  factors  at  the  end  of  the  increment. 

3.6  Fatigue  Analysis 

In  the  early  1960’s,  it  was  demonstrated  that  fracture  mechanics  is  a  useful  tool 
for  investigating  fatigue  crack  propagation.  In  this  section,  we  discuss  the  fatigue 
analysis  of  a  laser- welded  joint  under  constant  amplitude  cyclic  loading.  Consider  a  crack 
under  constant  amplitude  cyclic  load  AM.  Based  on  the  assumption  of  stress  intensity 
factor  being  directly  proportional  to  the  applied  load,  the  stress  intensity  amplitude  AK  at 
the  crack  tip  has  a  cyclic  variation  as  depicted  in  Fig.  3.5.  If  the  cyclic  plastic  zone 
formed  at  the  crack  tip  is  small,  the  crack  growth  rate  can  be  defined  by  the  expression: 

£  =  fW,  R).  (26) 

where  AK  =  (Kmax  -  Kmin)  is  the  amplitude,  R  =  Kmin/Kmax  is  the  load  ratio,  e  is  the  crack 
extension  measured  from  the  notch  tip  in  the  undamaged  configuration  and  de/dN  is  the 
crack  growth  per  cycle.  A  typical  fatigue  crack  growth  curve  for  metals  is  illustrated  in 
Fig.  3.6  on  a  log-log  plot  of  de/dN  versus  AK.  The  sigmoidal  curve  has  three  regions. 
The  high  and  low  AK  levels  show  non-linearity,  while  in  the  intermediate  region  or 
region-II,  the  curve  is  linear.  At  the  low  end,  when  AK  value  is  less  than  the  threshold 
value  AK,h,  de/dN  approaches  zero,  which  means  that  the  crack  will  not  grow.  For  high 
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Figure  3.5:  Constant  amplitude  stress  intensity  factor  at  the  crack  tip  due  to 

applied  cyclic  load. 


Figure  3.6:  Typical  fatigue  crack  growth  behavior  in  metals. 
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AK  values,  the  most  widely  accepted  hypothesis  is  that  crack  growth  rate  accelerates  as 
Kmax  approaches  the  fracture  toughness  of  the  material  K/c.  For  intermediate  values  of 
AK,  namely  within  region-II,  the  crack  growth  rate  can  be  described  by  the  power  law 
expression  of  Paris  and  Erdogan  (1960), 

%  =  C  *Km.  (27) 

where  AK  is  amplitude  of  the  cyclic  stress  intensity  factor,  and  C  and  m  are  material 
constants  which  are  obtained  experimentally.  Studies  have  shown  that  the  exponent  m 
ranges  from  2.0  to  7.0  for  most  materials.  Equation  (27)  is  widely  referred  to  as  the  Paris 
Law.  After  Paris’s  work  a  number  of  researchers  have  developed  equations  to  model  all 
or  part  of  the  sigmoidal  curve.  Some  researchers  proposed  a  relationship  for  region-II 
and  III,  assuming  that  region-III  is  caused  by  superposition  of  fracture  and  fatigue  rather 
than  plastic  zone  effects.  A  number  of  equations  attempt  to  describe  crack  growth  in  all 
three  regions.  A  simplified  form  of  one  such  equation  given  by  Forman  (1967)  is: 

de  CAKm 

dN  ~  ((1  -  R)KC  -  AK)  (28) 

At  region-II,  equation  (28)  will  reduce  to  equation  (27).  Generally  these 
equations  contain  four  material  constants:  C,  m,  Kc,  and  AK,h.  None  of  these  equations 
have  history  dependence,  and  thus,  are  only  valid  for  constant  (stress  intensity)  amplitude 
loading. 

The  objective  of  the  fatigue  crack  propagation  analysis  is  to  predict  the  fatigue 
life  of  the  joint.  This  can  be  defined  as  the  number  of  cycles  required  for  the  specimen  to 
fail  at  prescribed  stress  amplitude.  FRANC2D  finite  element  program  is  used  to 
propagate  the  crack  and  obtain  the  history  of  stress  intensity  factors  at  each  crack  tip.  All 
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crack  growth  analyses  are  carried  out  until  a  crack  extends  very  close  to  a  boundary  or 
another  crack  surface. 


A  MATLAB  program  is  written  to  compute  the  crack  extension  as  a  function  of 
the  number  of  cycles  using  the  stress  intensity  factor  histories  obtained  from  FRANC2D. 
This  is  achieved  by  integrating  the  Paris  Law.  The  power  law  exponent  m  in  the  Paris 
law  is  chosen  to  be  that  of  mild  steel  m  =  3.0  (Barsom  1971).  The  non-dimensional 
results  presented  in  the  following  sections  are  independent  of  the  material  constant  C  in 
the  Paris  law.  The  mathematical  steps  involved  in  this  formulation  are  summarized 
below.  Consider  the  propagation  of  the  crack  that  originates  from  the  right  notch  of  a 
joint  that  is  subjected  to  a  positive  cyclic  moment.  The  following  expressions  for  the 
mode  I  and  II  stress  intensity  factor  were  presented  previously, 


(29) 


where  the  subscript  i  signifies  the  mode.  For  any  arbitrary  moment,  AK  depends  on  AM, 
h  and  the  crack  extension.  It  should  be  noted  that  for  fatigue  analysis,  the  non- 


dimensional  stress  intensity  factor  /iM+,r  is  also  a  function  of  the  crack  extension  e/h. 


which  is  different  from  non-welded  overlap  length  a/h.  For  fatigue  crack  propagation,  a 
particular  joint  geometry  needs  to  be  focused  on  and  thus  other  dimensions  (w/h,  H/h  and 
a/h)  are  considered  to  be  fixed.  The  above  equation  can  thus  be  simplified  as, 


(30) 


Substitution  of  equation  (30)  into  equation  (27)  gives 


(31) 
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which  when  integrated  yields 


1  /7t3/2\m  fe°  de 

C\m)  'Jo  /(0 


dN  =  /V, 


(32) 


Equation  (34)  can  be  expressed  in  a  more  simplified  form  as, 

i  *(¥«)  r»  «© 

w  C '  (AM)m  '  J0  fQm'  (33) 

where,  zlA/  is  the  moment  amplitude.  The  above  equation  (equation  (33)),  is  then 
numerically  integrated  using  an  adaptive  simpson  quadrature  technique.  The  factor  / 
obtained  at  each  step  by  applying  a  unit  moment  AM  to  a  joint  of  unit  thickness  h,  is 
really  the  stress  intensity  factor  AK,  i.e.,  AK  =  f(e/h)  from  equation  (30).  The  stress 
intensities  obtained  at  each  crack  increment  has  been  recorded  and  while  considering 
mixed  mode  fracture  conditions,  effective  stress  intensity  factor  ( AKeg )  is  calculated  from 
them.  The  following  effective  stress  intensity  factor  is  used  since  it  gives  good 
correlation  with  experimental  results  (Qian  and  Fatemi,  1996), 

A  Keff  =  [A  K*  +  8.  A/f„4]1/4,  (34) 

The  total  crack  extension  range  is  divided  into  n  segments,  and  at  each  segment  AKeg  is 
evaluated  by  performing  cubic  interpolation.  The  AKeff  is  normalized  with  respect  to  the 
applied  load  and  upper  plate  thickness  ( h )  as  given  among  the  normalized  equations  in 
Section  3.4. 

Six  joint  configurations  are  selected  to  be  studied  for  fatigue  life  predictions,  of 
which,  one  is  set  to  be  common  among  each  loading  case.  Plots  of  the  normalized 
number  of  cycles  versus  normalized  crack  extension  are  shown  for  each  loading  case. 
There  normalized  number  of  cycles  for  the  moment  and  tensile  loads  are  defined  as, 
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N  =  (AM)mN, 


(35) 


where  AM  = 


C1/mAM 

3  i  t  and 


h'2  mj 


N  =  (A  P)mN, 


(36) 


—  C '  Hr  _ 

where  AP  =  — j-t- •  By  normalizing  the  number  of  cycles,  the  N  —  e  plot  becomes 

h(2+mJ 

independent  of  the  material  constant  C,  load  AM  or  AP,  and  upper  plate  thickness  h. 

In  the  manual  crack  growth  propagation  analysis  performed  here,  stress  intensities 
at  both  the  notches  are  obtained  and  an  effective  stress  intensity  factor  is  evaluated 
(equation  (34))  at  every  crack  increment.  One  tip  at  a  time  is  selected  and  an  increment 
size  for  each  crack  tip  is  calculated  using  the  following  procedure  and  is  supplied  to  the 
program.  Let  AKi  and  AK2  denote  the  effective  stress  intensity  factors  at  the  1st  and  2nd 
crack  tips  of  the  joint,  respectively.  For  a  specified  increment  of  load  cycles  AN  the  crack 
extensions  Aei  and  Ae2  will  be  different.  It  is  safe  to  assume  here  that  the  SIF  of  one 
crack  tip  would  be  greater  than  the  other,  say  AKt  >  AK2.  The  crack  with  the  largest 
stress  intensity  factor  is  allowed  to  propagate  by  the  maximum  crack  increment  Aei  =  Ae. 
Each  step  of  the  crack  propagation  analysis  corresponds  to  a  fixed  number  of  cycles,  i.e., 
AiVj  =  AN2  =  AN.  On  the  basis  of  Paris  Law, 


(37) 


(38) 


Dividing  equation  (37)  by  equation  (38), 


(39) 
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This  explains  that  the  extension  of  crack  2  is  a  fraction  of  the  extension  of  crack  1. 
Comparing  this  with  equation  (25),  it  can  be  inferred  that  power  b  is  actually  the  material 
constant  m  of  the  Paris  Law.  Therefore,  a  power  value  of  b  =  3.0  is  used  so  that  the 
results  obtained  from  the  crack  propagation  analyses  can  be  further  utilized  for  fatigue 
life  prediction. 

3.7  Results  and  Discussion 

Three  loading  cases  are  presented  in  this  section,  namely,  positive  moment, 
negative  moment  and  tensile  force.  Two  joint  configurations  are  discussed  in  each 
loading  case.  Among  the  two  joints  of  each  case,  one  is  set  to  be  common  to  ease  the 
purpose  of  comparing  the  performance  of  one  joint  under  different  loading  conditions. 

3.7.1  Loading  Case-I  (positive  moment  M*) 

To  ascertain  the  dependence  of  the  crack  propagation  trajectory  on  the  maximum 
crack  growth  size  Ae,  three  analyses  with  different  normalized  values  of  Ae  =  0. 1 ,  0.05 
and  0.025  were  performed  under  positive  moment.  These  analyses  were  run  just  to  get  an 
idea  of  the  increments  which  should  be  used.  The  analyses  are  run  to  a  maximum 
normalized  crack  extension  of  e  =  0.9  in  each  case.  By  analyzing  the  results  based  on  the 
crack  trajectory,  convergence  of  stress  intensity  factor  values,  and  the  computation  time 
taken  by  each  of  the  analysis,  a  standard  crack  increment  of  1 0%  of  the  weld-width  ratio 
w/h  is  selected  for  further  study. 
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3.7.1. 1  Common  Joint  (JC:1-0.6-1) 

Manual  crack  propagation  analysis  is  carried  out  on  the  joint  JC:1 -0.6-1  with  a 
unit  positive  cyclic  moment.  The  trajectory  obtained  is  illustrated  in  Fig.  3.7.  It  is 
observed  that  initially  both  the  left  and  right  cracks  begin  to  propagate.  After  few  steps 
the  left  crack  does  not  propagate  further,  since  the  computed  crack  growth  at  each  step  is 
smaller  than  the  chosen  tolerance  of  0.002. 

A  comparative  plot  of  the  normalized  effective  stress  intensity  factors  at  the  right 
and  left  notches  versus  the  corresponding  crack  extensions  e  is  shown  in  Fig.  3.8. 
Effective  stress  intensity  factor  values  are  plotted  to  emphasize  the  fact  that  the  joint  fails 
under  mixed  mode  condition.  The  dipping  down  of  the  left  notch  curve  shows  that  the 
left  crack  ceases  to  propagate  further  due  to  reduction  in  the  effective  stress  intensity 
factor  value. 

The  failure  load  for  each  crack  tip  is  computed  using  the  maximum  tensile  stress 
fracture  criterion  during  the  crack  propagation  process  and  the  minimum  of  the  two  is 
chosen  as  the  residual  fracture  strength  of  the  entire  joint.  The  notation  Mc  can  thus  be 
called  as  the  residual  strength  of  the  joint  at  a  certain  step  during  the  crack  propagation 
process,  since  it  is  the  maximum  load  the  joint  can  bear  just  before  crack  propagation  or 
material  failure.  Depicted  in  Fig.  3.9  is  a  plot  of  the  normalized  residual  strength  as  a 
function  of  the  normalized  crack  extensions  of  each  notch.  The  plot  is  very  useful  for 
evaluating  the  residual  strength  of  the  joint  as  the  cracks  propagate  through  the  structure. 
For  a  given  value  of  the  normalized  crack  extension  of  the  right  crack  tip,  the 
corresponding  value  of  the  residual  strength  can  be  inferred  from  Fig.  3.9.  For  example, 
the  residual  strength  and  crack  extension  of  the  left  notch  corresponding  to  a  right  notch 
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Figure  3.7:  Crack  propagation  trajectory  for  M+  load  for  the  joint  geometry  (H/h 

w/h  =  0.6,  a/h  =  1). 


e 

Figure  3.8:  Normalized  effective  stress  intensity  factor  vs.  normalized  crack 
extension  under  M  +  load  at  both  the  left  and  right  notches. 
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^ _  Residual  strength 
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Figure  3.9:  Normalized  residual  strength  of  the  joint  vs.  normalized  crack 

extension  under  M  +  load. 


Figure  3.10:  Normalized  N-e  curves  for  the  left  and  right  tips  under  a  positive 

moment  amplitude. 
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crack  extension  of  e=0.4  are  Mc=  0.151  and  e,=0.04,  respectively.  It  is  observed  that  the 
residual  strength  of  the  joint  increases  a  bit  in  the  initial  stages  of  crack  growth  and 
reaches  a  peak  value  after  which  it  decreases  monotonically.  This  is  because  the  mode  II 
stress  intensity  factor  is  significant  in  the  undamaged  joint.  In  the  initial  steps  of  crack 
propagation,  the  mode  II  stress  intensity  factor  decreases  significantly  while  the  mode  I 
stress  intensity  factor  increases  by  only  a  small  amount.  This  results  in  a  decrease  in 
effective  stress  intensity  factor  and  therefore  an  increase  in  the  residual  strength. 
Subsequently,  the  mode  I  stress  intensity  factor  continues  to  increase  thus  causing  the 
residual  strength  to  decrease  monotonically. 

The  normalized  crack  extensions  of  both  notches  versus  the  normalized  number 
of  cycles  obtained  by  integrating  the  Paris  law  are  plotted  in  Fig.  3.10.  The  left  notch 
propagates  only  up  to  e  =  0.04.  In  comparison  the  right  notch  keeps  growing  until  it 
reaches  close  to  the  top  surface  of  the  joint  where  e  =  0.87.  The  crack  growth  accelerates 
rapidly  as  it  approaches  the  upper  boundary  due  to  the  sharp  increase  in  stress  intensity 
factor  AK.  The  right  crack  reaches  close  to  the  upper  surface  at  N  =  0.0088  cycles.  One 
of  the  benefits  of  plotting  the  normalized  quantities  is  that  it  becomes  possible  to  obtain 
the  physical  crack  extension  e  if  the  number  of  cycles  N,  cyclic  moment  amplitude  AM, 
joint  thickness  h  and  material  constant  C  are  specified.  For  example,  if  the  applied  cyclic 
moment  is  sufficiently  small,  the  crack  will  extend  to  the  upper  surface  when  the 
normalized  number  of  cycles  is  0.0088.  The  normalization  equation  (35)  can  be  used  to 
determine  the  corresponding  number  of  cycles  N  as  follows, 


N  =  0.0088 


3m 

h  2 


+1 


C.  (AM)m‘ 


114 


For  example,  if  A M=  1  kip-in/in,  m  =  3,  h  =  0.5”  and  C  =  3.6x1  O'10  (in  IPS  units 
with  AK  specified  in  ksi-in  and  e  in  inches  in  the  Paris  Law)  the  right  crack  will  extend 
close  to  the  upper  surface  (e  =  0.435”)  at  N=  0.54  million  cycles.  This  will  be  elaborated 
upon  in  the  joint  case  studies  at  the  end  of  the  Chapter.  Another  benefit  of  presenting  the 
normalized  results  is  that  a  single  curve  is  sufficient  to  compute  the  fatigue  failure  limit 
of  a  joint  at  different  load  amplitudes  (AM)  and  upper  plate  thicknesses  (h).  Therefore, 
the  most  desired  curve  AM  -  N  (also  called  as  S-N  curve)  for  one  type  of  joint,  can  be 
obtained  with  the  help  of  just  one  plot.  The  normalized  moment  amplitude  has  already 
been  defined  in  equation  (35). 

If  the  cyclic  moment  magnitude  is  small,  failure  of  the  joint  can  be  defined  as  the 
point  at  which  the  crack  extends  close  to  a  boundary  or  another  crack  surface.  In  this 
case,  the  normalized  number  of  cycles  N ,  corresponding  to  maximum  crack  extension 
can  be  used  in  equation  (35)  to  obtain  the  normalized  moment  AM  as  a  function  of  the 
actual  number  of  cycles  N.  The  resulting  AM  vs  N  curve  is  a  straight  line  when  plotted 
on  the  log-log  axes  as  shown  in  Fig.  3. 1 1 . 

3.7.1. 2  Joint  Configuration  (JC:1  -0.2-1) 

Another  joint  configuration  JC:l-0.2-l  is  also  studied  due  to  its  susceptibility  of 
crack  growth  at  both  the  notches  under  positive  moment  load.  The  joint  geometry  is  the 
same  as  the  common  joint  studied  earlier,  except  for  the  weld  width,  which  has  been 
reduced  to  w/h  =  0.2.  The  details  of  the  joint  and  the  corresponding  crack  trajectory  are 
shown  in  Fig.  3.12.  Initially,  the  right  crack  grows  with  maximum  increments  while  the 
left  one  propagated  with  smaller  ones.  Due  to  propagation  of  right  crack  towards  the 
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Figure  3.11:  Log-log  AAf-N  (or  fatigue  life)  plot  with  positive  moment  amplitude 

for  the  joint  JC:l-0.6-l. 


Figure  3.12:  Crack  propagation  trajectory  for  M  +  load  for  the  joint  ( H/h  =  1, 

w/h  =  0.2,  a/h  =  1). 
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other  crack  tip  through  the  weld,  within  a  few  steps  the  rate  of  propagation  of  the  left 
crack  also  increases  with  an  increase  in  the  increment  size.  The  analysis  is  terminated 
when  the  left  crack  tip  approached  the  right  crack  surface.  It  is  observed  that  JC:1 -0.2-1 
has  a  very  localized  failure  in  the  welded  region. 

A  comparative  plot  of  normalized  effective  stress  intensity  factors  at  the  right  and 
left  notches  versus  e  is  shown  in  Fig.  3.13.  As  explained  above,  the  plot  illustrates  that 
initially  the  right  crack  had  higher  stress  concentration  and  that  it  grows  with  maximum 
supplied  increments.  Although  the  left  notch  propagated  with  a  smaller  crack  increment 
calculated  on  the  basis  of  equation  (39),  it  shows  a  higher  growth  rate  which  is  the  effect 
of  high  stress  concentration  developed  due  to  two  very  close  crack  tips. 

The  curve  shown  in  Fig.  3.14  depicts  the  normalized  residual  strength  of  the  joint 
at  certain  normalized  crack  extensions.  It  is  observed  that  the  residual  strength  of  the  joint 
decreases  rapidly  in  the  initial  stages  of  crack  growth. 

Fig.  3.15(a)  illustrates  the  normalized  number  of  cycles  versus  normalized  crack 
extension  curves  for  both  the  notches  in  the  joint  considered.  The  convergence  of  the  red 
curve  at  N  =  0.0034  signifies  that,  at  this  N,  the  left  crack  grows  steeply  and  meets  all 
the  way  to  the  right  crack  surface  which  results  in  the  failure  of  the  joint. 

As  done  before,  the  normalized  number  of  cycles  to  failure  N  =  0.0034  is  used  in 
conjunction  with  equation  (35)  to  obtain  the  fatigue  life  curve  shown  in  Fig.  3.15(b).  A 
comparison  of  N  at  maximum  crack  extension  for  joints  JC:1 -0.6-1  and  JC:1 -0.2-1 
shows  that  the  fatigue  life  N  for  a  given  cyclic  load  decreases  by  61.36%  if  the  weld- 
width  is  decreased  from  w/h=0.6  to  wlh-02. 
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Figure  3.13:  Normalized  effective  stress  intensity  factor  vs.  normalized  crack 
extension  under  M '+  load  at  both  the  left  and  right  notches. 


extension,  under  positive  moment. 
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(a) 


Figure  3.15:  (a)  Normalized  N-e  curves  for  the  left  and  right  tips,  and  (b)  Fatigue 
life  plot,  for  JC:l-0.2-l,  under  positive  moment  amplitude. 
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3.7.2  Loading  case-ll  (Negative  Moment,  li/T) 

Crack  propagation  analysis  is  also  carried  out  for  a  negative  moment  load  case. 
Two  joint  configurations  JC:l-0.6-l  and  JC:2.5-0.2-l  are  selected  to  be  analyzed  under 
this  loading  condition. 

3.7.2.1  Common  Joint  (JC:1-0.6-1) 

On  the  common  joint  geometry  (JC:l-0.6-l),  16  manual  propagation  steps,  with  a 
power  value  of  3,  are  performed.  Fig.  3.16(a)  illustrates  the  crack  trajectory  at  both  the 
right  and  left  notches.  The  right  crack  trajectory  shows  a  small  dip  down,  which  signifies 
the  presence  of  hinging  effect  discussed  in  Section  2.6.5.  However,  it  has  been  observed 
that  for  lower  H/h  ratios  the  hinging  effect  is  not  very  prominent.  This  is  because  a  lower 
plate  of  smaller  thickness  will  be  less  stiff  and  thus  will  provide  less  reaction  force. 

Variations  in  the  normalized  effective  stress  intensity  and  the  normalized  residual 
strength  with  respect  to  the  crack  extension  e  have  been  shown  in  Fig.  3.16(b)  and  (c). 
The  crack  trajectory  and  the  stress  intensity  factor  plots  are  quite  different  from  that  of 
the  positive  moment  loading  case  (see  Figs.  3.7  and  3.8).  In  this  case  the  opening  is  on 
the  opposite  side  (left),  and  the  crack  trajectories  at  both  the  tips  are  downwards  due  to 
the  effect  of  negative  moment.  Careful  inspection  of  the  right  crack  trajectory  shows  that 
it  initially  propagates  upward  and  then  changes  direction  and  extends  downwards.  This  is 
due  to  the  stresses  caused  by  the  hinging  effect  in  the  initial  stages.  Due  to  the  change  in 
the  propagation  angle  of  the  right  crack  tip,  the  residual  strength  is  affected  and  thus  it 
shows  some  variation  in  the  initial  steps.  A  little  bump  can  be  seen  on  the  effective  stress 
intensity  factor  at  the  left  tip  in  Fig.  3.16(a)  in  the  initial  stages.  The  increase  in  the 
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(a) 


(c) 


intensity  factor  and  (c)  Normalized  residual  strength  vs.  normalized  crack 
extension  e,  under  a  M'  load  in  the  joint  geometry  JC:l-0.6-l. 
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effective  stress  intensity  factor  at  the  left  tip  during  the  initial  steps  is  due  to  the  changes 
in  the  mode-II  stress  intensity  factor.  Since  after  a  certain  crack  extension,  the 
propagation  becomes  more  and  more  mode-I,  the  Kn  effects  recede  and  the  curve 
smoothens. 

Figure  3.17(a)  shows  the  left  and  right  crack  extensions  as  a  function  of  the 
normalized  number  of  cycles.  It  is  seen  that  the  left  crack  rapidly  grows  and  reaches  the 
lower  boundary  when  N  =  0.002.  The  corresponding  fatigue  life  curve  is  given  in  Fig. 
3.17(b).  A  comparison  of  N  at  maximum  crack  extension  from  Figs.  3.10  and  3.17(a) 
shows  that  the  direction  of  the  moment  load  has  a  significant  effect  on  the  fatigue  life  of 
the  joint.  The  fatigue  life  N  of  the  joint  under  a  negative  moment  is  77.3%  smaller  than 
when  it  is  subjected  to  a  positive  moment  of  identical  cyclic  moment  amplitude. 

The  converged  normalized  number  of  cycles  obtained  from  Fig.  3.17(a),  is  again 
utilized  to  predict  the  fatigue  life  of  the  joint.  For  this  purpose,  the  normalized  negative 
moment  amplitude  AM  applied  on  the  joint  is  plotted  against  actual  number  of  cycles  N  to 
failure  of  the  joint,  on  a  log  scale.  Fig.  3.17(b)  simply  explains  that  if  a  certain  amount  of 
AM  is  applied  on  the  joint,  then  the  joint  will  fail  in  the  corresponding  number  of  cycles 
based  on  the  plot. 
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(a) 


Normalized  Number  of  Cycles  N  x  10 


(b) 


AM-N  plot,  for  negative  moment  amplitude,  for  the  joint  JC:l-0.6-l. 
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3.7. 2.2  Joint  Configuration  (JC:2.5-0.2-1) 

A  detailed  joint  model  with  dimensions  H/h  =  2.5,  w/h  =  0.2  and  a/h  =  1  is  made 
and  meshed  in  CASCA.  Fig.  3.18(a)  illustrates  the  crack  trajectory  for  this  model.  It  can 
be  seen  that  the  crack  propagates  through  the  weld,  presenting  a  fairly  complicated  crack 
propagation  trajectory.  One  has  to  be  careful  when  choosing  the  maximum  crack 
increment  Ae  during  the  crack  propagation  analysis.  If  the  increment  is  too  small,  the 
analysis  will  require  an  excessive  number  of  elements  while  if  the  increment  is  too  large 
the  trajectories  are  very  coarse  and  it  leads  to  the  premature  intersection  of  the  two  crack 
boundaries. 

Normalized  values  of  effective  stress  intensity  factor  history  and  residual  strength 
with  respect  to  normalized  crack  extension  are  plotted  for  each  step  in  Fig.  3.18(b)  and 
(c).  Fig.  3.18(b)  shows  that  the  effective  stress  intensity  factor  is  higher  at  the  left  tip 
during  the  initial  steps. 

The  normalized  residual  strength  at  e  =  0  in  this  case,  is  found  a  little  higher  than 
the  previous  joint  case  JC:l-0.6-l.  By  comparing  Fig.  3.18(c)  and  Fig.  3.16(c),  one  can 
say  that  at  a  fixed  a/h  =  1,  increasing  H/h  produces  a  joint  of  higher  strength  although  the 
weld- width  w/h  is  smaller.  Conversely,  the  smaller  the  weld  length  of  the  joint,  the  more 
significant  drop  is  seen  in  its  strength  with  crack  growth. 

For  the  joints  under  negative  moment  amplitude,  e.g.  JC:  1  -0.6-1,  some  amount  of 
hinging  effect  has  been  seen.  This  effect  becomes  more  prominent  with  either  an  increase 
in  the  non-welded  overlap  length  or  decrease  in  the  welded  length.  This  is  explained  with 
a  reference  of  another  joint  JC:2.5-0.2-l  and  its  N  —  e  plot  shown  in  Fig.  3.19(a).  The 
curves  shown  in  Fig.  3.19(a)  and  Fig.  3.17(a),  clearly  shows  the  effect  of  hinging 
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(a) 


(c) 


intensity  factor  and  (c)  Normalized  residual  strength  vs.  normalized  crack 
extension  e ,  under  a  M '  load  in  the  joint  geometry  JC:2.5-0.2-l. 
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(a) 


AM-N  plot,  for  negative  moment  amplitude,  for  JC:2.5-0.2-l. 
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mechanism  in  both  the  negative  moment  cases.  In  Fig  3.17(a),  the  left  notch  propagated 
at  a  faster  rate  due  to  comparatively  smaller  non-welded  and  higher  welded  length.  While 
in  Fig.  3.19(a),  right  notch  showed  a  higher  growth  rate  due  to  high  stress  concentration 
as  a  result  of  hinging. 

Fig.  3.19(b)  shows  the  fatigue  life  prediction  of  joint  JC:2.5-0.2-l.  The  end  point 
of  the  line  is  defined  as  fatigue  endurance  limit  of  the  joint.  This  means  that  the  crack  will 
not  grow  at  any  applied  load  below  the  lowest  point  of  the  line.  By  comparing  the  fatigue 
life  plots  of  the  two  joint  in  this  loading  case,  one  can  say  that  by  increasing  the  H/h  ratio 
the  fatigue  endurance  limit  of  the  joint  increased.  However,  it  did  not  increase 
significantly  due  to  the  reduced  welded  length.  Therefore,  it  is  worth  adding  here  that 
JC:2.5-0.6-l  under  negative  moment  amplitude  has  70%  higher  fatigue  life  than 
JC:1 -0.6-1. 

3.7.3  Loading  Case-Ill  (Positive  Tension,  P *) 

Referring  to  the  future  goal  of  the  study,  if  the  bolted  and  riveted  joints  in  a  ship 
structure  are  to  be  replaced  with  laser  welded  joints,  then,  it  would  be  desirable  to 
understand  the  behavior  of  the  joint  under  a  tensile  load.  Keeping  this  in  mind,  a  study 
has  been  conducted  to  investigate  crack  growth  under  a  tensile  load. 

3.7.3.1  Common  Joint  (JC:1-0.6-1) 

On  the  common  joint-geometry  (JC:1  -0.6-1)  a  constant  distributed  axial  load  is 
applied  at  the  overhang  edge  of  the  upper  plate.  With  a  power  value  of  3,  and  16  steps  of 
propagation,  the  crack  trajectory  obtained  is  almost  the  same  as  shown  in  Fig.  3.16(a). 
However,  they  are  not  identical.  This  is  because  in  this  case,  no  crack  propagation  is  seen 
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on  the  right  notch  (Fig.  3.20(a)),  which  explains  that  hinging  effect  is  really  seen  only  in 
the  negative  moment  cases. 

The  normalized  A^and  normalized  Mc  plots  shown  in  Fig.  3.20(b)  and  (c)  are  for 
left  notch  of  the  joint  JC:  1  -0.6-1.  Corresponding  normalized  N-e  curve  and  fatigue  life 
plots  are  shown  in  Fig.  3.21(a)  and  (b).  These  plots  are  self  explanatory. 
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100 


(C) 


intensity  factor  and  (c)  Normalized  residual  strength  vs.  normalized  crack 
extension  e  at  the  left  notch,  under  a  P  +  load  in  the  joint  geometry  JC:l-0.6-l. 
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(a) 


Figure  3.21:  (a)  Normalized  N-e  curves  for  the  left  tip,  and  (b)  Log-log  AM-N 
plot,  under  tensile  force  amplitude,  for  JC:l-0.6-l. 
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3.7. 3. 2  Joint  Configuration  (JC:2.5-0.2-1) 

Another  joint  type  is  considered  for  this  loading  case  in  order  to  check  the  effect 
of  varying  w/h  and  H/h  on,  the  crack  trajectory,  the  stress  intensity  factor  and  the  critical 
load,  for  an  applied  tensile  load.  Fig.  3.22(a),  show  the  crack  trajectory  for  the  joint 
JC:2.5-0.2-l.  It  is  observed  that,  although,  effective  stress  intensity  factor  values  at  both 
the  notches  in  this  case  are  quite  close,  the  cup  and  cap  trajectory  obtained  here  is  really 
the  effect  of  decrease  and  then  increase  in  the  mode-II  stress  intensity  factor,  which 
affected  the  crack  propagation  angle  too.  It  can  be  seen  that  for  a  smaller  w/h  ratio,  even 
if  the  lower  plate  is  kept  thick,  the  crack  propagates  through  the  weld  when  the  specimen 
is  subjected  to  a  tensile  load.  The  variation  in  the  normalized  effective  stress  intensity 
factor  and  the  normalized  residual  strength  for  this  joint  with  respect  to  the  crack 
extension,  have  been  shown  in  Figs.  3.2(b)  and  (c).  As  plotted  earlier,  the  red  color  lines 
depict  the  left  and  blue  depicts  the  right  notches  for  the  joint  geometry.  It  is  interesting  to 
note  that  there  comes  a  propagation  step  when  both  the  crack  tips  have  nearly  the  same 
effective  stress  intensities.  At  that  step  the  crack  increments  calculated  for  both  the  crack 
tips  are  almost  the  same. 

Fig.  3.23(a)  and  (b)  shows  the  fatigue  plots  of  normalized  number  of  cycles 
versus  normalized  crack  extension,  and  normalized  applied  tensile  force  amplitude  versus 
actual  number  of  cycles  to  complete  failure  of  the  joint,  respectively.  There  is  not  much 
difference  in  the  fatigue  life  of  the  two  joint  types  under  this  loading  case,  but  the  one 
with  higher  H/h  ratio  can  bear  more  number  of  cycles  at  the  same  applied  load. 
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22 


e 


(c) 


intensity  factor  and  (c)  Normalized  residual  strength  vs.  normalized  crack 
extension  e,  under  a  P +  load  in  the  joint  geometry  JC:2.5-0.2-l. 
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Normalized  Applied  Moment  A P  Normalized  Crack  Extension  e 


(a) 


Normalized  Number  of  Cycles  N  x  10 


Figure  3.23:  (a)  Normalized  N-e  curves  for  the  left  and  right  tips,  and  (b)  Log-log 
AAf-N  plot,  for  tensile  force  amplitude. 
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3.8  Joint  Case  Studies 


An  effort  has  been  done  in  the  previous  sections  to  elaborate  on  the  variation  of 
stress  intensity  factors,  residual  strength  and  number  of  cycles  with  crack  extension,  but 
the  results  are  presented  in  a  non-dimensional  form  that  can  be  applied  to  a  variety  of 
joint  configurations.  Two  case  studies  are  performed  in  this  section  to  elucidate  on  the 
implementation  of  those  results  and  get  some  data  for  practical  applications.  Designers 
are  interested  in  the  residual  strength  and/or  fatigue  life  of  a  structure.  In  continuation,  it 
is  also  desirable  to  know  the  residual  strength  of  the  joint  after  it  has  been  subjected  to  a 
certain  number  of  cycles  under  a  prescribed  cyclic  load  magnitude. 

The  first  case  study  is  performed  on  the  joint  model  JC:1 -0.6-1  with  a  repeated 
cyclic  load  {A\f)  where  Mmm  =  0.  It  is  assumed  that  the  joint  is  made  of  steel  with  C  = 
3.6  x  10'10,  m  =  3,  Kic=  109  ksiVin,  AAT,h  =  5  ksiVin  and  it  has  an  upper  plate  thickness  of 
h  =  0.25”.  The  Paris  law  constants  C  and  m  used  in  the  following  case  studies  correspond 
to  that  of  a  mild  steel  (Barsom  1971).  The  output  of  the  codes  used  for  fatigue  plots,  is 
used  to  generate  a  plot  of  residual  strength  versus  normalized  number  of  cycles.  Fig.  3.24 
contains  a  plot  of  the  residual  strength  of  the  joint  versus  the  normalized  number  of 
cycles.  The  plot  points  out  the  minimum  and  maximum  residual  strength  of  the  joint 
0.269  and  Mcmax=  4.639  kips-in/in)  at  zero  and  maximum  crack  extensions  at  one 
of  the  notches.  It  is  interesting  to  see  that  the  plot  shows  some  hardening  effects  up  to 
certain  N  which  corresponds  to  a  certain  crack  extension,  after  which  the  strength  of  the 
joint  drops  almost  linearly.  The  little  peak  in  the  residual  strength  plot  experienced  here, 
has  been  referenced  back  to  the  crack  propagation  analysis  section  of  this  chapter,  and  it 
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Figure  3.24:  Normalized  number  of  cycles  versus  strength  left  in  the  joint 
JC:l-0.6-l  for  an  applied  positive  moment  amplitude. 


Figure  3.25:  Applied  moment  amplitude  vs.  number  of  cycles  for  JC:l-0.6-l, 
considering  failure  of  the  structure. 
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is  observed  that  this  peak  is  a  result  of  the  large  mode-II  stress  intensity  factor  at  the  very 
first  step  of  crack  propagation. 

A  plot  of  moment  amplitude  versus  number  of  cycles  AM-N  is  considered  as  a 
major  result  of  any  fatigue  life  analysis.  It  has  been  experimented  and  mentioned  many 
times  in  the  literature  that  AM-N  plot  for  steel  shows  an  ‘L’  shape,  if  broader  range  of 
applied  moments  and  number  of  cycles  to  failure  are  considered.  This  essentially  means 
that  there  is  a  fatigue  endurance  limit  of  the  structure  below  which  it  will  not  fail  for  any 
number  of  cycles.  An  attempt  has  been  made  to  relate  this  plot  to  the  sigmoidal  curve 
discussed  above.  For  a  better  representation  of  Region-II  of  the  sigmoidal  curve  shown  in 
Fig.  3.6,  it  is  required  to  define  the  lower  and  upper  limits  of  the  applied  load  amplitude. 
Critical  moment  or  residual  strength  calculations  have  already  been  done  to  get  the  upper 
limit  of  the  Region-II  of  the  sigmoidal  curve.  The  threshold  load  AM,hreshoid-  0.2128  kips- 
in/in  is  evaluated  using  Maximum  Tensile  Stress  Criterion  and  AKth  of  the  material. 
Several  moment  amplitudes  are  applied  on  the  model  between  AM  threshold  and  some  30% 
below  AMcmax,  and  number  of  cycles  to  failure  of  the  joint  is  recorded,  as  shown  in  the 
Fig.  3.25.  For  fatigue  analysis  applied  loads  are  not  carried  up  to  the  critical  value,  but  are 
kept  below  that  value.  The  curve  mimics  the  S-N  curve  for  any  other  arbitrary  specimen 
and  is  meant  for  the  failure  of  the  structure  and  not  for  the  failure  of  the  material. 

The  second  design  study  is  performed  on  JC:2.5-0.2-l  joint  with  tensile  load  P  \ 
An  upper  plate  thickness  h  =  0.25”  is  used  and  same  formulation  and  code  is  used,  except 
using  hm  :  instead  of  h3m'2  for  normalized  quantities  in  the  moment  cases.  The  output 
obtained  by  running  the  code  are  illustrated  in  Fig.  3.26(a)  and  (b). 
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(a) 


Figure  3.26:  (a)  Normalized  number  of  cycles  vs.  strength  left ,  for  an  applied 
tensile  force  amplitude  (b)  Applied  tensile  force  amplitude  vs.  number  of  cycles  to 

failure  of  the  joint  JC:2.5-0.2-l. 
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3.9  Conclusions 


In  this  chapter,  the  finite  element  code  FRANC2D  is  used  to  analyze  the  crack 
growth  trajectories  of  different  joint  configurations  and  to  compute  the  fatigue  life  by 
performing  fatigue  analysis  of  those  joints  under  different  loads.  It  is  observed  that  for 
different  joint  configurations  under  different  types  of  loading  conditions,  the  crack 
trajectories  are  significantly  different.  Some  situations  are  harder  to  explain,  especially 
where  there  are  crack  growths  at  both  the  notch  tips.  In  such  conditions,  one  crack  starts 
growing  towards  another  which  significantly  affects  the  stress  intensity  around  the  other 
crack.  These  different  trajectories  also  explain  that  the  same  joint  behaves  differently 
under  different  loading  conditions,  and  therefore  a  crack  trajectory  of  any  stake 
connection  cannot  be  predicted  for  any  arbitrary  load.  However,  one  can  say  that  the 
crack  would  always  tend  to  grow  through  the  weld  when  the  weld-width  is  small. 

The  manual  crack  propagation  gives  better  control  over  the  operation  of  the 
program.  In  performing  a  manual  propagation,  one  can  keep  the  power  value  constant 
(which  according  to  the  multiple-crack  fatigue  crack  growth  method  is  really  the  material 
constant  m )  and  compute  the  increments  at  each  crack  tip.  The  crack  incrementation  has 
been  kept  above  a  tolerance  of  0.002,  which  means  the  crack  is  not  propagated  if  the 
increment  calculation  gives  below  0.002  for  a  crack  tip.  For  the  geometries  like  the  one 
studied  here,  mixed  mode  stress  intensity  factors  are  required  to  be  considered. 
Therefore,  referring  to  the  power  law  crack  growth  model,  the  stress  intensity  factor 
amplitude  AK  is  considered  as  a  mixed  mode  parameter  AKeff.  The  above  computed 
increments  are  then  supplied  to  the  program  to  get  new  crack  tips.  The  analysis  is  then 
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run  with  dynamic  relaxation  algorithm  for  a  tolerance  of  lxl O'4,  and  stress  intensities  are 
recomputed  at  the  new  crack  tips. 

The  method  discussed  above,  to  generate  the  normalized  plots,  would  turn  out  to 
be  very  handy  when  a  designer  would  have  to  compute  the  strength,  number  of  cycles  it 
can  withstand  and  the  fatigue  life,  of  a  joint,  which  already  has  some  crack  growth.  For 
example,  if  a  structure  with  n  cracks  is  considered  and  the  total  crack  extension  at  any  of 
the  crack  tips  in  that  structure  is  measured,  then  with  the  help  oiMc  —  e  and  N  —  e  plots, 
it  will  be  easy  to  compute  the  residual  strength  of  the  structure  and  the  maximum  number 
of  cycles  it  can  withstand. 

For  fatigue  analysis  six  joint  configurations  are  analyzed  using  FRANC2D.  A 
simple  form  of  fatigue  evaluation  is  done  without  experiments,  all  based  on  numerical 
methods  by  using  Paris  Law.  In  an  attempt  to  visualize  the  performance  of  the  joints 
under  a  broad  range  of  applied  loads,  maximum  and  threshold  values  are  also  considered. 
By  assuming  a  threshold  value  of  the  stress  intensity  factor  (AK,h  =  5%  of  AKC),  applied 
threshold  load  amplitude  ( AM,h )  is  computed.  These  maximum  and  threshold  values  are 
utilized  to  predict  the  fatigue  life  of  each  joint  configuration  considered.  It  also  helped  in 
explaining  that  the  crack  will  not  grow  at  low  loads,  and  thus  defines  the  fatigue 
endurance  limit  of  a  particular  joint  under  consideration.  Two  joint  case  studies  are 
performed  to  clarify  the  normalization  concept  carried  over  throughout  the  chapter. 
These  case  studies  develop  a  better  understanding  of  converting  the  normalized  results 
into  practical  use. 

It  has  been  seen  in  the  literature  (Anderson,  1995)  that  for  more  ductile  materials, 
stress  intensity  approach  would  not  be  too  appropriate  and  that  AK  in  the  Paris  Law 
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needs  to  be  replaced  by  AJ  (a  J-integral  amplitude  load).  To  add  to  this,  before  applying 
AJ,  one  may  realize  that  the  joint  geometries  discussed  in  this  study  produces  mixed¬ 
mode  crack  propagation  conditions.  This  area  has  a  great  potential  for  current  and  future 
research. 
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